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INTRODUCTION 


Although  the  finite  element  method  is  a proven  effective  method 
for  obtaining  numerical  solutions  of  solid  mechanics  problems1, 
its  impact  on  computational  fluid  dynamics  has  been  felt  only  in  the 
past  few  years2.  Because  of  the  diverse  applications  of  this  method  in 
continuum  mechanics,  many  departures  from  the  original  method  used  in 
structural  analysis  have  been  made.  Our  adaptation  of  the  finite 
element  method  for  direct  application  to  unsteady  gas  flows  in  two 
spatial  variables  is  based,  in  part,  on  Lynn  and  Arya's  least  squares 
formulation3’1*  and  on  Polk's  one  dimensional  study5.  Lynn  and  Arya's 
approach  is  based  on  the  elementwise  least  squares  minimization  of 
the  differential  residual  error  which  allows  a direct  finite  element 
formulation  from  the  governing  differential  equations.  Furthermore, 
since  the  governing  equations  are  hyperbolic,  the  finite  elements 
can  be  constructed  in  both  space  and  time  so  that  they  approximate  the 
domain  of  determinancy  associated  with  hyperbolic  problems.  Polk  com- 
bined these  two  concepts  and  applied  them  to  the  unsteady  isentropic 
flow  of  an  inviscid  gas  expanding  behind  a piston.  Using  both  linear 
and  quadratic  approximations  to  the  dependent  variables,  he  showed 
good  agreement  between  the  finite  element  results  and  the  exact  solu- 
tion for  the  smooth  portion  of  the  flow.  Near  the  gradient  discontin- 
uities which  occurred  in  the  flow,  Polk  constructed  the  finite  elements 
so  that  a side  of  the  element  and  the  locus  of  discontinuities  coincided. 
Using  this  special  construction,  good  agreement  was  obtained  everywhere 
in  the  flow. 

This  report  .is  concerned  with  a finite  element  method  for  unsteady 
inviscid  compressible  flows  in  two  spatial  dimensions,  a corresponding 
pilot  computer  code  and  some  resulting  numerical  experiments. 


^ O.C . Zienkiewicz , The  Finite  Element  Method  in  Engineering  Science, 
McGraw-Hill 3 1971. 

2 

J.T.  Oden , O.C.  Zienkiewicz , R.H.  Gallagher , C.  Taylor , eds .,  Finite 
Element  Methods  in  Flow  Problems , University  of  Alabama  in  Huntsville 
Press,  1974. 

SP.P.  Lynn  and  S.K.  Ary a,  "Use  of  the  Least  Squares  Criterion  in  the 
Finite  Element  Formulation ,"  Int.  J.  Hum.  Meth,  Engng. , 6,  75-88 , 1972. 

4 

P.P.  Lynn  and  S.K.  Arya,  "Finite  Elements  Formulated  by  the  Weighted 
Discrete  Least  Squares  Method , " Int.  J.  Hum.  Meth.  Enging.,  8,  71-90 , 
1974. 

5J.  Polk,  "A  Least  Squares  Finite  Element  Approach  to  Unsteady  Gas 
Dynamics,"  BEL  Report  No.  1885,  May  1976.  (AD  #A026531) 
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The  system  of  governing  equations  are  nonlinear  hyperbolic  equa- 
tions. We  take  advantage  of  this  fact  to  simplify  the  general  finite 
element  methodology.  This  is  contrary  to  the  parabolic  regularization 
method  of  Oden,  et.  al.6  for  hyperbolic  problems.  In  the  techniques  of 
Oden,  et  al.,  certain  terms  which  depend  on  the  discretization  para- 
meters are  appended  to  the  equations  so  that  they  become  parabolic. 

This  parabolic  problem  is  then  solved  by  a finite  element  technique. 

It  can  be  shown  for  a class  of  problems  that  the  solution  of  the 
parabolic  problem  converges  to  the  original  hyperbolic  solution  in  the 
limit  as  the  mesh  size  tends  to  zero.  On  the  other  hand,  our  formu- 
lation deals  directly  with  the  hyperbolic  equations. 

Our  construction  of  the  finite  elements  is  reminiscent  of  Polk's 
construction  in  that  they  are  in  both  space  and  time  but  differ  in  that 
they  enclose,  not  coincide  with,  the  domain  of  dependence.  It  will  be 
shown  that  this  construction  simplifies  the  necessary  integration 
routines  while  satisfying  a Courant  condition.  No  special  construction 
of  the  finite  element  is  provided  near  steep  gradients  and  discontin- 
uities in  the  flow.  Consequently,  no  special  knowledge  of  the  solution 
is  required  and  all  interior  nodes  in  the  calculation  are  treated 
identically.  Furthermore,  by  applying  Lynn  and  Arya's  least  squares 
minimization  to  each  finite  element,  we  can  avoid  the  large  matrices 
generally  associated  with  the  finite  element  method  for  elliptic  prob- 
lems while  still  retaining  the  essential  advantages  of  the  method. 

The  general  methodology;  the  construction  of  the  finite  elements, 
the  approximations  to  the  flow  variables  and  the  formulation  in  terms 
of  the  least  squares  error  criterion,  is  explained  in  Section  II. 

Section  III  contains  a brief  discussion  of  the  computer  code,  the 
form  of  the  governing  equations  and  certain  approximations  used  within 
the  code.  The  results  of  two  numerical  experiments  are  given  and 
discussed  in  Section  IV.  Section  V contains  a summary  of  the  method  and 
areas  in  which  future  work  is  required. 


II.  GENERAL  METHODOLOGY 

The  first  step  in  the  finite  element  methodology  is  to  divide  the 
solution  region  into  elements.  We  divide  the  computational  domain  for 
a given  time  (a  two  dimensional  region)  into  triangular  elements.  The 
vertices  of  the  triangles  are  called  nodes.  We  assume  that  the  bound- 
aries are  stationary  so  that  the  triangular  divisions  remain  unchanged 


J.T.  Oden,  L.C.  Well ford,  and  C.T.  Reddy,  "A  Study  of  Convergence  and 
Stability  of  Finite  Element  Approximation  of  Shock  and  Acceleration 
Waves  in  Nonlinear  Materials,"  U.S.  Army  Research  Office  Report  P-11860- 
M/DAAG29-76-G-C022,  August  1976. 
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with  time.  The  system  of  governing  equations  for  compressible  fluid 
flow  include  coupled  nonlinear  partial  differential  equations  which  ex- 
press conservation  of  mass,  momentum,  and  energy  plus  an  algebraic 
equation  of  state.  Because_the  governing  equations  are  hyperbolic, 
the  solution  at  a point  (x,y,t+At)  in  the  solution  domain  depends  only 
on  the  value  of  the  flow  variables  within  the  intersection  of  the  domain 
of  dependence  (the  mach  cone)  from  the  point  (x,y,t+At)  and  the  (x,y,t) 
plane.  Consequently,  given  a_  unison  of  triangles  in  the  (x,y,t)  plane 
with  a vertex  at  the  point  (x,y,t)  and  the  values  of  the  flow  variables 
at  the  corresponding  nodes^  we_can  compute  a value  of  At  such  that  the 
mach  cone  from  the  point  (x,y,t+At)  lies  within  the  union  of  triangles. 
Since  the  computed  value  of  At  will  vary  from  node  to  node,  we  take  the 
minimum  value  of  At  over  all  nodes  as  the  next  time  step.  This  value 
allows  a systematical  advance  in  time. 

We  now  define  our  finite  element  at  a point_(x,^,t+At)  as  the  union 
of  all  prisms  with  a base  vertex  at  the  point  (x,y,t)  and  with  a uni- 
form height  At.  See  Figure  1.  The  present  one  offers  several  desirable 
simplifications  over  other  possible  constructions  of  the  finite  element. 
From  the  discussion  above  it  is  clear  that  the  values  of  the  flow  vari- 
ables at  a node  are  independent  of  the  values  at  the  other  nodes  at  the 
same  time  level.  Thus,  the  interconnection  of  the  nodal  values  which 
is  characteristic  of  finite  element  formulations  of  elliptic  problems 
and  which  result  in  the  manipulation  of  large  matrices  can  and  will  be 
avoided.  We  will  solve  for  the  central  nodal  values  at  the  new  time 
level  by  considering  only  the  finite  element  at  the  central  node. 
Futhermore,  any  necessary  time  integrations  over  the  finite  element  are 
simplified,  since  the  sides  of  the  elements  are  independent  of  time. 

The  next  step  is  to  choose  the  interpolating  or  trial  function  over 
the  finite  element.  Let  w be  a flow  variable;  that  is,  a dependent 
variable  computed  directly  from  the  governing  differential  equations, 
not  the  equation  of  state.  The  interpolating  function  for  w is 

w(x,y,t)  = uo(x,y,t)  + t*(a.x  + a.+1y  + ai+2) , (1) 

where  the  points  (x,y,t)  lie  within  the  finite  element  at  (x,y,t+At)  and 
the  parameters  a^,  a^  . , a.+2  are  to  be  determined.  The  function 
w (x,y,t)  represents1tne  flow  variable  w at  time  t and  is  assumed  known. 
THese  interpolation  functions  form  over  each  finite  element  a linear 
approximation  in  space  to  the  time  derivative  of  u>. 

To  complete  the  model,  a technique  to  determine  the  parameters  a. 
and  thus  the  flow  variables,  is  needed.  A basic  idea  in  the  finite  1 
element  methodology  is  to  minimize  the  errors  arising  from  the  residual 
of  the  governing  differential  equations  in  terms  of  the  interpolating 
functions.  Thus,  the  finite  element  method  approximates  the  minimum 
residual  whereas  the  finite  difference  method  approximates  the  differ- 
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ential  equation.  Since  the  proposed  analysis  is  elementwise,  we  desire 
a minimization  technique  for  each  element.  To  this  end,  we  choose 
the  elementwise  least  squares  minimization  of  the  differential  residual 
error  employed  by  Lynn  and  Arya. 

We  substitute  the  interpolating  functions  into  the  k—  governing 
differential  equation,  make  the^result  dimensionless  and  denote  the 
resulting  residual  by  D^Cx.y.tja) , where  a is  the  vector  of  unknown 
parameters  a^.  Each  residual  is  dimensionless  so  that  no  individual 
Dt  will  numerically  dominate  the  least  squares  sum  of  all  the  residuals 
because  of  dimensional  disparities.  The  D^'s  are  explicit  algebraic 
functions  of  the  independent  variables  x,y,t  and  the  parameters  a.  The 
total  error  in  the  sense^of  least  squares  over  a particular  finite 
element  is  denoted  by  E(a)  and  is  given  by 

, , NOEQ 

E(a)  = J J J J.  Dk(x,y,t;J)  dxdydt,  (2) 

V 


where  NOEQ  is  the  number  of  governing  equations  and  V is  the  volume  of 
the  finite  element.  We  wish  to  minimize  E(a)  with  respect  to  the  a.'s. 
A necessary  condition  for  the  existence  of  a minimum  is  8E/8a^  = 0 <$r 


u u, 

dxdydt  = 0 , for  each  i (3) 
i 


Note  that  the  derivatives  SD^/Sa^  can  be  explicitly  calculated.  By 
solving  the  nonlinear  algebraic  system  of  equations  (3)  for  the  un- 
knowns a.,  we  can  determine  the  values  of  the  flow  variables  at  the 
nodal  point  (x,"y,"t+At) . We  repeat  this  process  for  each  interior 
node  in  the  solution  domain. 

For  a boundary  node,  the  above  procedure  is  slightly  altered.  We 
rewrite  the  given  boundary  condition(s)  at  the  boundary  node  in  terms 
of  the  interpolating  functions  (1).  We  then  solve  for  the  unknown 
parameters  a.  in  equation  (3)  at  the  boundary  node  subject  to  the  re- 
written boundary  conditions.  Thus,  at  a boundary  node  we  no  longer 
have  a pure  minimization  problem  as  at  an  interior  node,  but  rather  a 
constrained  minimization  problem. 


The  pilot  computer  code  is  written  in  a modular  structure  fashion  in 
order  to  clarify  the  logic  of  the  program  and  to  allow  changes  in  the 
form  of  the  governing  equations,  the  interpolating  functions,  the  method 
selected  to  solve  the  nonlinear  system,  etc..  Such  flexibility  is  highly 
desirable  for  this  pilot  code.  For  example,  in  the  finite  difference 
techniques,  the  formulation  of  the  equations  have  a profound  affect  on  a 
method's  performance  (see,  for  example,  Moretti 7) . Similarly,  the  form 
of  the  equations  may  affect  the  performance  of  the  finite  element  method. 
Different  forms  of  the  equations  for  unsteady  compressible  flows  are 
listed  in  Appendix  A. 

The  code  has  three  major  components.  The  first  component,  sub- 
routine START,  accepts  the  geometric  and  control  parameters.  Further- 
more, this  section  accepts  and/or  generates  the  nodal  positions  within 
the  x-y  solution  subdomain  and  necessary  initial  and  boundary  value 
data.  The  second  component,  subroutine  TIMEf rEP,  calculates  the  time 
increment  (subroutine  DELTAT)  and  the  flow  variables'  values  at  each 
node  at  the  new  time  (subroutine  CALCI  for  r terior  nodes  and  CALCB 
for  boundary  nodes).  The  "heart"  of  the  pi*ot  code  is  clearly  the 
TIMESTEP  routine,  since  the  general  method  outli.  ed  in  Section  II  is 
implemented  in  this  portion  of  the  code.  The  third  component,  sub- 
routine DISPLAY,  provides  the  output  and  graphics  capabilities  for  the 
program.  The  code  as  listed  in  Appendix  D uses  the  non-conservation 
inviscid  form  of  the  governing  equations  in  Cartesian  coordinates, 
where  body  forces,  heat  absorption  and  heat  fluxes  are  neglected  (see 
equation  system  (Al)  of  Appendix  A),  the  Newton-Raphson  iteration 
method  and  certain  approximations.  We  briefly  discuss  this  specific 
situation  below. 

The  governing  equations  in  dimensional  variables  are: 


3p 

at 

8p 

+ u ax 

9p 

+ V ~“ 
9y 

+ P 

(4) 

p( 

3u  ^ 

at  u 

3u 

_—  + v 
3x 

3u  \ 
ay  / 

•U-o. 

(5) 

p( 

3v 

at  + u 

3v 

+ V 

3x 

av  \ 

ay  / 

. |£=  0, 
iy 

(6) 

( 

3e 

3e  ^ 

3e\ 

/au  av\ 

(7) 

p ( 

at  + u 

t—  + V 
3x 

3y)+ 

? ia3T  + a?)  = °* 

P 

= P(p,e)  . 

(8) 

Gino  Moretti,  "A  Pragmatical  Analysis  of  Discretization  Procedures  for 
Initial  - and  Boundary  - Value  Problems  in  Gas  Dynamics  and  Their  Influ- 
ence on  Accuracy  or  Look  Ma,  No  Wiggles! , " Polytechnic  Institute  of 
New  York  Report  No.  74-15,  September  1974. 
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Here  the  spatial  coordinates  are  x and  y,  t is  the  time,  u is  the 
x-component  of  the  velocity,  v is  the  y-component  of  the  velocity, 
p is  the  pressure,  p is  the  density  and  e is  the  internal  energy  per 
unit  mass.  The  functional  form  of  the  equation  of  state  is  given  by 
equation  (8).  The  particular  equation  of  state  used  in  a given  cal- 
culation is  specified  in  the  subroutine  EQNST. 

In  th£  actual  calculations  for  a given  finite  element  at  the 
point  (x,y,t+At),  we  translate  the  origin  of  the  coordinate  system 
to  the  point  (x,y,t)  (see  subroutine  TRANS L) . Conceptually,  the 
translation  enables  each  finite  element  to  have  its  base  centered  at 
the  same  point  (x,y,t)  = (0,0,0)  and  practically,  it  simplifies  the 
calculations.  The  interpolating  functions  for  the  variables  u,  v,  p, 
e within  a finite  element  are: 


u(x,y,t)  = uQ(x,y,0)  + t • (a^  + a2y  + a3) , (9) 
v(x,y,t)  = vo(x,y,0)  + t • (a4x  + a5y  + a&) , (10) 
P(x,y,t)  = po(x,y,0)  + t • (ayx  ♦ agy  + ag) , (11) 
e(x,y,t)  = eQ(x,y,0)  + t • (a1Qx  + any  + a12) , (12) 


where  the  variables  x,  y,  t are  in  the  translated  finite  element,  the 
variables  with  subscript  zero  denote  the  variables  at  the  known  time 
level  t=0,  and  the  a^'s,  i=l , 2,  ...,  12  are  the  unknown  parameters. 

We  define  the  dimensionless  residuals  Dj<.(x,y,t;'S) , k=l,2,3,4  as 
the  quantities  obtained  by  substituting  equations  (8) -(12)  into  the 
four  differential  equations  (4) -(7),  respectively,  and  then  by  dividing 
the  first  result  by  a ratio  of  a characteristic  density  to  At,  the 
second  and  third  results  by  a ratio  of  a characteristic  velocity  to  At, 
and  the  fourth  by  a ratio  of  a characteristic  internal  energy  per  mass 
to  At.  The  values  of  the  characteristic  quantities  are  the  density, 
sound  speed  and  internal  energy  per  unit  mass  at  the  center  node  and 
At=0.  The  derivatives  of  p with  respect  to  the  variables  x,  y and  t 
can  be  obtained  by  the  chain  rule.  The  subsequent  partial  derivatives 
of  p with  respect  to  p and  e are  calculated  directly  from  the  equation 
of  state  and  are  coded  into  the  subroutine  EQNST.  The  functions  u , 
vq,  p , e are  estimated  by  a linear  approximation  on  each  triangle 
based°on  ?heir  values  at  the  vertices.  Although  other  approximations 
are  possible,  the  present  one  is  easily  computable  and  is  independent 
of  the  number  of  prisms  in  the  finite  element.  The  partial  derivatives 
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of  uq,  vq,  pQ,  e , required  in  D^(x,y,t;t)  vary  from  triangle  to 
triangle°and°con§equently  equation  (3)  must  be  rewritten  as 


FjCa) 


3E(a)  WRl 

8ai  " di 


= 0 


III 

Jlth  prism 
i=l ,2, . . .12, 


4 

I 

k=l 


8Dk 

aT  dxdydt 

i 


(13) 


where  NOEQ  is  now  four  and  the  finite  element  contains  NUMTRI  prisms 
(or  triangles) . 


The  Newton-Raphson  method  is  used  to  solve  the  system  of  nonlinear 
algebraic^equations  (13)  (see  subroutine  SILVER)  ^ince  the  terms 
Dk(x»y»t;a)  ate  known  functions  of  the  parameter  a.  The  second  partial 
derivatives  of  the  least  squares  error  can  be  explicitly  calculated  and 
are  given  by 


tth  prism 

, i,j=l,2,. . . ,12,  (14) 


(see  subroutine  FU) . The  initial  values  for  a , a , a , a.  in  the 
Newton-Raphson  scheme  for  the  finite  element  at  (0,0,  At)  are  taken  as 
the  previously  computed  values  of  these  parameters  for  the  finite 
element  at  (0,0,0),  (see  subroutine  INITZR) . However,  at  the  first 
time  step  a special  calculation  is  needed  to  find  the  appropriate 
initial  values  (see  subroutine  INITNG) . These  parameters  are  determined 
by  solving  for  3/3t  in  each  of  the  equations  (4) -(7),  by  approximately 
the  spatial  dependence  of  the  variables  at  the  initial  time  by  linear 
functions  on  each  triangle,  by  finding  the  corresponding  values  of 
3/3t  at  each  vertex,  by  equating  these  to  the  time  partials  of  equations 
(9)  - (12),  by  averaging  the  resulting  values  of  the  a^'s  over  the  tri- 
angles and  by  translating  the  result.  In  all  cases,  the  initial  values 
of  a^,  a2,  a4,  ag,  a^,  ag,  a1Q,  and  a^  aTe  set  to  zero.  In  the  sample 
problems,  the  iteration  converged  faster  with  the  zero  values  than 
with  the  more  obvious  choice  a^  j*  0,  ij*3,6,9,  12.  The  integrals 
of  Dfc"  (B^Dj^/Sa^Saj ) , (3D^/3aj)  and  (3D^/3a^)  • (3Dj^/3aj ) are 

approximated  by  a two  point  Gaussian  quadrature  in  time  and  by  a pro- 
duct of  first  order  functions  in  space.  For  example. 
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where  t , m = 1,2,  are  the  two  Gaussian  quadrature  points  between  zero 
and  At.  In  the  spatial  integrals  at  each  Gaussian  time,  both  factors 
are  approximated  by  first  order  functions  in  x and  y which  are  then 
multiplied  and  integrated  exactly  over  the  desired  triangle  (see  sub- 
routines DSUB,  EVALPR  and  AREAIN) . 

Once  the  values  of  the  parameters  a.  are  known,  the  desired  values 
of  the  flow  variables  at  the  center  node  at  the  new  time  (t  = At)  can 
be  determined  easily  from  equations  (9)  - (12),  (see  subroutine  VALNEW) . 

The  above  discussion  applies  equally  to  interior  and  boundary  type 
nodes.  However,  as  noted  at  the  end  of  Section  II,  the  minimization  at 
a boundary  node  is  a constrained  minimization.  The  type  of  constraints 
will  depend  on  the  type  of  boundary  condition  imposed  at  a given  node. 

As  an  example,  the  zero  normal  velocity  boundary  condition  imposed  at  a 
solid  wall  is  incorporated  into  the  equation  system  (13)  in  Appendix  B. 


IV.  NUMERICAL  EXPERIMENTS 


In  this  section  the  results  of  two  non-steady  flow  calculations, 
the  flow  behind  a cylindrical  blast  wave  and  the  flow  across  a prop- 
agating normal  shock,  are  presented.  These  two  examples  are  computed 
in  order  to  ascertain  the  scheme's  characteristics  on  interior  nodes 
for  a smooth  and  discontinuous  flow,  respectively.  Both  of  these  time- 
dependent  flows  are  essentially  one  dimensional  in  space,  however,  they 
will  be  treated  as  two  dimensional  problems.  Furthermore,  since  closed- 
form  solutions  are  known  for  each  problem,  the  accuracy  of  the  finite 
element  formulation  can  be  precisely  evaluated. 


The  similarity  solution  of  the  blast  wave  problem  is  discussed  by 
Sedov0  and  the  formulas  for  the  cylindrical  case  are  summarized  in  Appen- 


L.I.  Sedov , Similarity  and  Dimensional  Methods  in  Mechanics,  Academic 
Press,  1959.  ' " * ” ' ’ 
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dix  C.  The  flow  behind  a cylindrical  blast  wave  rather  than  a planar 
wave  is  computed  because  of  its  associated  circular  solution  domain. 

The  computational  domain  for  a given  time  consists  of  an  annulus  from 
radius  ra  to  radius  rb  which  is  the  position  of  the  shock  front  at  the 
initial  time  t0.  Perhaps  the  most  important  advantage  of  the  finite 
element  method  is  its  capability  of  dealing  with  complex  geometrical 
shapes  by  using  arbitrarily  shaped  simple  elements.  The  ease  by  which 
this  cylindrical  problem  is  handled  by  the  Cartesian  finite  element 
program,  especially  the  circular  boundaries,  demonstrates  this  advan- 
tage in  an  elementary  manner.  In  Figure  2 the  first  quadrant  of  a 
computational  domain  for  a given  time  is  drawn,  where  ra  = 2.2[m], 
rb  = 3.0[m],  the  nodes  are  equally  spaced  at  four  degree  intervals  on 
a given  circle  and  the  radial  divisions  are  computed  so  that  approxi- 
mate equilateral  triangles  result.  Consequently,  the  size  of  the  tri- 
angles increases  with  the  radial  distance  from  the  origin.  Note  the 
good  approximation  to  the  arc  boundaries  by  the  series  of  straight 
lines  forming  the  base  of  the  triangles.  A finite  difference  method 
in  Cartesian  coordinates  either  would  approximate  the  arc  boundaries 
by  a series  of  horizontal  and  vertical  lines  or  would  require  a trans- 
form of  the  solution  domain.  In  both  cases  special  treatment  of  the 
boundaries  would  be  required.  On  the  other  hand,  no  special  program- 
ming is  needed  to  treat  the  arc  in  the  finite  element  code. 

Numerical  calculations  were  performed  for  the  case  ra  = 2.2[m], 
rb  = 3.0[m]  and  the  initial  time  t = 3.85342034xl0~3[s] . On  the 
initial  time  plane  the  values  of  flow  variables  are  calculated  from  the 
exact  solution  (see  Appendix  C) . For  the  inflow  condition  at  ra  and 
outflow  condition  at  rb  we  again  assign  the  exact  value  to  the  flow 
variables.  The  computed  ratios  of  p/ps,  vr/vTs,  e/es  and  p/ps  (sub- 
scripts denotes  value  at  the  shock  front)  are  compared  to  Sedov's 
exact  values  (solid  line)  in  Figures  3,  4,  5,  and  6,  respectively.  The 
symbols  1X1  and  ^ denote  the  computed  value  on  the  triangular  finite 
element  mesh  with  nodes  spaced  at  two  degree  intervals  (the  computed 
radial  subdivisions  range  from  0.07[m]  to  0.09[m])  and  at  four  degree 
intervals  (the  computed  radial  subdivisions  range  from  0.14[m]  to 
0.18[m]),  respectively.  Both  sets  of  values  are  at  the  same  time 
(4 . 12516608xl0"3  [s] ) . The  former  used  four  timesteps  (At  35  0.68xl0"^[s]) 
and  the  latter  two  timesteps  (At5®*  1 .36x10"^  [s] ) to  reach  the  termination 
time.  The  maximum  absolute  value  of  the  percent  relative  error  for  the 
ratios  of  pressure,  radial  velocity,  internal  energy  and  density  are 
1.15%,  0.12%,  0.10%  and  1.24%,  respectively,  for  the  finer  mesh  and 
2.68%,  0.30%,  0.36%  and  3.04%,  respectively,  for  the  coarser  mesh.  We 
recall  that  the  density  is  computed  from  the  equation  of  state  once  the 
pressure  and  energy  are  calculated  from  the  differential  equations. 

Hence,  the  error  in  the  density  includes  both  the  pressure  and  energy 
errors.  We  note  that  the  maximum  relative  error  in  each  flow  variable 
occurs  where  the  finite  elements  are  the  largest;  that  is,  near  rb.  At 
the  opposite  end,  near  ra,  the  finite  elements  are  the  smallest  and  the 
absolute  value  of  the  percent  relative  error  of  the  pressure,  radial 
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rs  = 3.10397969  (m) 
ps  = 1.48919895  * 105  (Pa) 


Figure  3.  Comparison  of  Exact  Normalized  Pressure  for  Cylindrical 
Blast  Wave  [8]  with  Computed  Values  for  Nodes  Spaced 
at  4°  (I)  and  2°  (n)  Intervals 


NORMALIZED  POSITION  (r/rs) 


r$  = 3.10397969  (m) 

vr$  =3.  13520549  *102  (m/$) 


Figure  4.  Comparison  of  Exact  Normalized  Radial  Velocity  for 
Cylindrical  Blast  Wave  [8]  with  Computed  Values  for 
Nodes  Spaced  at  4°  (X)  and  2 (m)  Intervals 
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Figure 


rs  = 3.10397969(m) 
e$  =4.9  1 475  6 74  xiO4  (J/kg) 


5 Comparison  of  Exact  Normalized  Internal  Energy  for 
Cylindrical  Blast  Wave  [8]  with  Computed  Values  for 
Nodes  Spaced  at  4°  (X)  and  2°  (m)  Intervals 
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rs  = 3.10397969  (m) 

Ps  = 7.57514068  (lc8/m3) 


Figure  6.  Comparison  of  Exact  Normalized  Density  for  Cylindrical 
Blast  Wave  [8]  with  Computed  Values  for  Nodes  Spaced 
at  4°  (X)  and  2°  (x)  Intervals 
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velocity,  internal  energy  and  density  ratios  for  the  finer  and  courser 
mesh  are  0.11%,  0.01%,  0.02%,  0.14%,  and  0.44%,  0.18%,  0.19%,  0.64%, 
respectively.  Hence,  overall  and  despite  the  relatively  large  size  of 
the  triangular  elements,  the  agreement  is  fairly  good. 

A normal  shock  of  strength  five  (p2/Pi  = 5)  propagating  into  a rec- 
tangular field  is  the  second  calculation.  The  subscripts  1 and  2 on  the 
variables  denote  their  value  in  the  pre-  and  post  shock  states,  respec- 
tively. The  flow  domain  consists  ideally  of  the  infinite  slab  of  the 
(x,y)  plane,  -°°  < x < °°,  Y2  .1  Y Yl  f°r  time  greater  than  the  initial 
time  tQ.  The  shock  is  initially  at  the  position  yS(),  y2  < ySo  < Yj  > 
and  is  moving  in  the  positive  y-direction.  The  quiescent  state  is  char- 
acterized by  zero  velocities  and  by  pressure  and  density  of  air  at  sea 
level  and  at  temperature  288. 15° [K] . The  values  of  the  flow  variables 
in  the  disturbed  region  are  computed  by  the  Rankine-Hugoniot  relations. 
The  formulas  and  values  used  to  calculate  the  exact  solution  are  listed 
in  Appendix  C.  The  nodes  are  equally  spaced  along  constant  y values 
and  the  y directed  divisions  are  computed  so  that  approximate  equilat- 
eral triangles  result.  The  computational  domain  for  a given  time  is 
very  similar  to  that  in  Figure  2 except  that  the  radial  and  angular 
variables  now  correspond  to  the  y and  x variables,  respectively.  How- 
ever, all  the  triangles  in  this  case  are  the  same  size.  On  the  initial 
time  plane  t0,  the  values  of  the  flow  variables  are  calculated  from  the 
exact  solution.  The  boundaries  at  y2  and  y^  are  taken  sufficiently  far 
away  from  the  shock  front  so  that  the  constant  values  in  both  the  dis- 
turbed and  quiescent  regions  are  obtained  by  the  flow  variables. 

For  the  calculations  the  y directed  subdivisions  are  small  in  order 
to  model  the  shock  with  its  extremely  thin  thickness.  The  y directed 
divisions  are  0.005[m],  y2  = 2.965[m],  yj  = 3.050[m],  ys  = 3.0025[m] 
and  t0  = 0.0[s].  The  graphs  of  the  computed  ratios  for  ?he  pressure 
(p/ps) , the  y-velocity  (v/vs) , the  internal  energy  (e/es)  and  the  density 
(p/ps)  versus  the  y-axis  are  given  in  Figures  7,  8,  9,  10,  respectively, 
for  three  times;  1 . 1524902x10"^ [s]  (fourth  time  level),  2.5986071xl0"4  [s] 
(ninth  time  level)  and  4. 0654310xl0"4 [s] (fourteenth  time  level).  The 
symbol  denotes  the  theoretical  position  of  the  shock  front.  The 

Figures  7,  8,  9,  and  10  show  fair  resolution  of  the  locations  of  the 
shock  front  for  the  given  times.  Spurious  oscillations  develop  at  the 
shock  front  as  expected,  since  no  artifice  was  introduced  into  the  equa- 
tions or  scheme  to  suppress  them.  The  oscillations  occur  about  the  exact 
solutions  p/ps  = v/vg  = e/eg  = p/p5  = 1 in  the  disturbed  flow  region. 

A major  objective  of  future  work  is  to  curtail  these  oscillations 
(see  Section  V of  this  report) . 

We  close  this  section  with  a discussion  of  the  calculation  time.  The 
times  needed  to  calculate  the  flow  variables  both  for  an  entire  time 
level  and  at  an  individual  interior  node  are  given  to  help  ascertain  the 
time  characteristics  of  the  scheme.  All  the  calculations  were  done  on 
the  BRLESC  computer  facility  at  the  Ballistic  Research  Laboratory.  The 
calculations  of  the  flow  behind  a cylindrical  blast  wave  took  approxi- 
mately 1.3  and  0.5  minutes  for  the  finer  and  coarser  grids,  respectively, 
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Figure  7.  Variation  of  Normalized  Pressure  with  Distance  for 
Normal  Shock  of  Strength  Five  at  Three  Times 
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t = 1.1524902  xi0'4(s) 


Figure  8.  Variation  of  Normalized  y Component  of  Velocity  with 
Distance  for  Normal  Shock  of  Strength  Five  at 
Three  Times 
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Figure  9. 


Variation  of  Normalized  Internal  Energy  with  Distance 
for  Normal  Shock  of  Strength  Five  at  Three  Times 
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Figure  10.  Variation  of  Normalized  Density  with  Distance  for 
Normal  Shock  of  Strength  Five  at  Three  Times 
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of  run  time  per  time  level  and  approximately  eight  seconds  per  interior 
node  for  both  grids.  These  times  include,  on  the  average,  two  iterations 
per  node  in  the  Newton-Raphson  method.  The  shock  propagation  calcula- 
tions took  approximately  4.5  minutes  per  time  level  and  approximately 
16  seconds  per  interior  node.  The  shock  propagation  problem  typically 
required  4 iterations  per  node  in  order  to  find  convergent  values  of 
the  unknown  parameters  a^.  Although  the  above  times  are  not  as  small  as 
one  might  hope,  they  can  be  reduced  substantially  by  simple  modifica- 
tions of  the  algorithm  (see  Section  V).  Such  simple  optimizations  must 
be  investigated  in  future  work. 


V.  SUMMARY 

We  have  presented  a numerical  scheme  for  solving  time  dependent 
hyperbolic  equations  in  two  space  dimensions.  This  method  merges  the 
concepts  of  the  finite  element  method  and  the  properties  of  hyper- 
bolic systems  of  equations  and  is  applied  to  unsteady  gas  flows.  The 
corresponding  hydrocode  is  summarized  and  listed.  Finally,  numerical 
calculations  involving  both  smooth  and  shocked  flows  are  given  and 
discussed . 

The  purpose  of  this  report  is  to  give  a summary  of  the  work  al- 
ready accomplished  rather  than  a definitive  description  of  the  quality 
and  usefulness  of  this  numerical  scheme.  However,  several  positive 
aspects  of  the  method  can  already  be  seen.  The  formulation  of  the 
method  is  straight-forward  and  avoids  the  large  matrices  associated 
with  the  finite  element  method.  The  method  handles  different  geomet- 
rically shaped  boundaries  easily  and  accurately.  Even  for  relatively 
large  mesh  sizes,  the  results  for  a smooth  flow  are  accurate.  Finally 
in  a propagating  shock  problem,  the  shock's  position  can  easily  be 
discerned . 

From  the  example  calculations  in  Section  IV,  it  is  clear  that 
improvements  must  be  made  in  several  areas  before  the  method's  po'-^n- 
tial  can  be  accurately  ascertained.  Listed  below,  not  necessar  i 
order  of  importance  or  difficulty,  are  several  such  areas. 

1.  Curtail  the  spurious  oscillations  near  shocks.  Several  meth- 
ods; such  as  the  flux-corrected  transport  techniques  of  Boris,  et  al,9 
and  the  well-known  artificial  viscosity  method  (see  Roache10),  can  be 


9 

D.L.  Book , J.P.  Boris , and  K.  Hain  "Flux-Control  Transport  II: 
Generalizations  of  the  Method ,"  J.  Comp.  Phys.,  18,  pp.  248-83, 

July  1975. 

10 

P.J.  Roache,  Computational  Fluid  Dynamics,  Hermosa  Publishers,  P.  0. 
Box  8172,  Albuquerque,  New  Mexico,  1976. 
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applied.  Although  the  latter  is  simpler  to  use,  the  former  technique 
may  hold  mere  promise  because  the  oscillations  do  occur  about  the 
correct  solution  and  the  flux  correction  will  not  significantly  reduce 
the  resolution  of  the  shock  as  does  artificial  viscosity. 

2.  Shorten  computing  time.  The  run  time  can  be  significantly  re- 
duced by  simple  alterations  in  the  algorithm.  Recall  that  the  spatial 
derivatives  of  the  variables  at  the  known  time  level  (u  , v , p0,  e ) 
were  computed  on  every  triangular  base  of  the  finite  elemen?  which  ° 
resulted  in  the  repeated  calculation  (up  to  NUMTRI  times)  of  the  resid- 
uals D]<  and  their  partial  derivatives  at  the  nodes  (see  equations  (13)  - 
(14)).  If  these  derivatives  at  the  nodes  were  computed  once  for  a given 
finite  element,  the  time  reduction  would  be  a factor  of  0.5.  The  run 
time  could  be  further  reduced  by  an  order  of  magnitude  if  the  calcula- 
tions were  done  on  a machine  similar  to  the  CDC  7600. 

3.  Rerun  examples  with  different  equation  formulations.  Certain 
formulations  may  increase  the  accuracy  of  the  computed  results  and 
decrease  the  oscillations  due  to  shocks.  See  Appendix  A for  different 
form  of  the  governing  equations. 

4.  Apply  the  method  to  actual  numerical  problem.  By  applying  this 
method  to  a problem  with  complicated  boundaries,  not  only  could  the 
method's  treatment  of  boundaries  be  tested,  but  also  the  entire  method 
could  be  compared  to  another  numerical  solution  technique. 

5.  Extend  the  method  with  respect  to  "infinite"  strength  shocks 
and  moving  boundaries.  Once  these  extensions  are  incorporated  into 
the  method,  this  scheme  could  hopefully  be  used  to  finally  develop  an 
adequate  model  of  the  severe  transitional  ballistics  environment. 
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APPENDIX  A 


VARIOUS  FORMULATIONS  OF  THE  GOVERNING  EQUATIONS 


The  system  of  governing  equations  for  fluid  flow  consists  of  dif- 
ferential equations  which  express  the  conservation  of  mass  (continuity 
equation),  momentum  (Navier-Stokes  equations)  and  internal  energy 
(energy  equation)  plus  an  algebraic  equation  of  state.  We  list  four 
different  formulations  of  these  time  dependent  differential  equations 
in  two  spatial  variables.  In  the  first  three  subsections  the  follow- 
ing definitions  are  used:  the  first  and  second  coefficients  of  viscos- 
ity, heat  absorption,  the  x and  y components  of  the  body  force  and  heat 
flux  are  denoted  by  y(x,y,t)  [kg/(m*s)],  A(x,y,t)  [kg/(m*s)],  H(x,y,t) 

[J/  (kg*s)  ] , B^x.y.t)  [N/kg] , B2(x,y,t)  [N/kg] , q^x.y.t)  [J/(m2-s)],  and 

q9(x,y,t) [J/(m^*s)] , respectively. 

1.  Independent  Cartesian  coordinates  are  x,y,t  and  dependent 
valuables  are  u,  v,  p,  p,  e. 


3p 

3t 


+ h (pu)  + h (pv)  = 0 


’I 


vavilbleTlTr!1*-1  £artesian  variables  are  x,y,t  and  dependent 
riables  are  u,  v_j_  e,  p,  p.  The  variables  O'  = pu,  v = pv  B,  = dB, 

!nH_hP52’/  = pe?  H = pH  are  the  x and  y components  of  the  momentum 
respectively?6 * intCrnal  ener^  and  heat  absorption  per  unit  volume, 
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3.  Independent  Cartesian  variables  are  x,y,t  and  dependent  vari- 
ables are  vr,  0,  S,  p,  p.  The  variables  vr (x,y,t) [m/s]  given  by  vr  = 

(u2  + v2)1/2  , 9(x,y,t) [rad]  given  by  tan  0 = v/u,  S[J/K],  and  T[K]  are 
the  modulus  of  the  velocity  vector,  the  argument  of  the  velocity  vector, 
entropy  and  temperature,  respectively. 


3p 

at 


+ COS0 


(PV  ‘ pVr  If 1 = ° 
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• n I 30  30  • a 3V 

-V  sxn0  — + V T—  sin©  - . 
r V3t  r 3y  3y 


[ Ij  <pvr>  * cvr  If]  * s1"0  [I? 

[f  3vr  3v  230  "N 

cos0  + vr  _I  cos0  . vr  ^ sin©  J 

— cos()^j 

(sir 

t (Vr  If  cosQ  + aT  sinG)  ] 

[(^T  - Vr  If)  sinG  + (ar  + Vr  i)  cos0] 
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4.  Independent  variables  are  r,z,t  and  dependent  variables  are 
Ur’  U0’  W’  P**  P*>  I*  In  the  governing  equations  listed  below,  the  fol- 
lowing definitions  apply:  r is  the  radial  distance  (r  ■ (x2  + y2)l/2) 

[m] , z is  the  axial  distance  [m] , t is  the  time  [s] , ur(r,z,t)  is  the 
radial  velocity  [m/s] , u (r,z,t)  is  the  swirl  velocity  [m/sj,  w(r,z,t) 
is  the  axial  velocity  [m/s] , p*(r,z,t)  is  the  density  [kg/mJ] , p*(r,z,t) 
is  the  pressure  [Pa]  and  I is  the  internal  energy  per  unit  mass  [J/kg] . 
Furthermore,  the  first  and  second  coefficients  of  viscosity,  the  radial 
and  axial  components  of  the  heat  flux  vector,  the  heat  absorption  term 
and  the  radial,  angular  and  axial  components  of  the  body  force  are 
denoted  by  v(r,z,t) [kg/(m-s)] , n (r , z ,t) [kg/ (m- s) ] , hr(r,z,t)(J/ (m2-s) ] , 

hz (r,z ,t) [ J / (m2*s) ] , H*(r , z ,t) [J/ (kg* s) ] , Br(r,z,t) [N/kg] , BQ(r,z,t) 

[N/kg],  Bz(r,z,t) [N/kg] , respectively.  In  this  formulation,  the  flow  is 
three  dimensional  and  axially  symmetric;  that  is,  the  six  unknowns  de- 
pend only  on  the  three  variables,  r,z,  and  t. 
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APPENDIX  B 


BOUNDARY  NODE  FORMULATION  FOR  A ZERO  NORMAL  VELOCITY  BOUNDARY  CONDITION 

Let  Q be  a node  at  a stationary  wall  where  the  tangent  is  defined. 
The  normal  velocity  being  zero  at  Q implies  that 

u(x,y,t)  sina  - v(x,y,t)  cosa  = 0 (Bl) 

where  (x,y)  are  the  spatial  coordinates  of  node  Q and  a is  the  angle 
between  the  tangent  line  at  Q and  the  x-axis.  Since  the  wall  is  sta- 
tionary a is  a constant  and_e£uation  (Bl)  holds  for  all  times^  By 
translating  the  axes  to  (x,y,t),  (the  plane  corresponding  to  t is  the 
initial  plane),  and  using  the  interpolating  functions  for  u and  v, 
equations  (9)  and  (10),  respectively,  equation  (Bl)  becomes 

[uq (0,0,0,)  + a^t]  sina  - [vq (0,0,0)  + a^t]  cosa  = 0,0<t<At.  (B2) 

Since  the  solid  wall  is  stationary,  equation  (B2)  reduces  to 

a-j  sina  - a^  cosa  = 0 . (B3) 

To  minimi ze_the_l east  square  residual  error  over  the  finite  element  at 
the  point  (x,y,t  + At)  subject  to  equation  (B3) , we  use  Lagrange  multi- 
pliers to  obtain 


cosa 


NUMTRI 

1 

1=1 


I Ini 


£th  prism 


dxdydt 


a^sina  - a^  cosa  = 0.  (B5) 

In  the  minimization  of  the  residual  error  over  the  boundary  type  ele- 
ment at  which  zero  normal  velocity  is  imposed,  equations  (B4)  and  (B5) 
replace  equations  dE(t)/da^  = 0 and  9E(i)/9a6  = 0 of  equation  system 
(13),  respectively. 


APPENDIX  C 


EXACT  SOLUTIONS  OF  TEST  PROBLEMS 

Cylindrical  Blast  Wave.  Sedov's  solution8  (pp.  219-20)  for  a 
cylindrical  blast  wave  generated  by  the  instantaneous  release  of  a 
finite  amount  of  energy  proportional  to  EQ  into  a gas  with  initial 
density  p0  is  summarized  below.  For  the  calculations,  the  specific 
heat  ratio  of  a perfect  gas  y and  two  constants  Eq  and  pQ  are: 

y = 1.40, 

E0  = 6.887025656E+06  [J/m] , (Cl) 

p0  = 1.262523446  [kg/m3]  . 

The  similarity  solution  which  uses  the  Rankine-Hugoniot  strong  shock 
conditions  can  be  expressed  in  terms  of  a parameter  8 as  follows: 

rs  = 48 .32791537*  (t)  1/2  [m]  , 

Vrs  973 . 1614183/rs  [m/s]  , 

ps  = 7.575140676  [kg/m3]  , 

ps  = 1.434797011xl06/rs2  [Pa]  , (C2) 

r/rs  = 0.9562806477(146  - S)1/?  (56  - 7$2)'l/2 , 

vr/vr  = 2.40  • 6 • r/r  , 

p/ps  = 2. 5 12993254x10* 4 • (146  - 5)5/?  [ (5-76)/ (1-26) ] 10/3, 

p/ps  = 6 . 61 8423368x10* 3 • 6 • [ (5-76)/ (1 -26) ] 7/3  , 

where  the  subscript  s denotes  the  value  at  the  shock  front  and  the  para- 
meter 6 satisfies  the  inequalities  5/14  £ 8 ^ 5/12.  Note  that  8 = 5/14 

corresponds  to  the  origin  of  the  blast,  a singularity  for  the  energy 
and  sound  speed  and  6 = 5/12  corresponds  to  the  position  of  the 
shock  front  at  a given  time.  The  internal  energy  for  a perfect  gas  of 
specific  heat  ratio  y = 1.4  is  given  by  e = 2.5  p/p.  The  x and  y compo- 
nents of  the  velocity  are  easily  computed  by  the  formulas  vrcos6  and 

v^sin  6,  respectively,  where  the  tan  6 = y/x. 


Propagating  Normal  Shock.  The  quiescent  state  into  which  the  nor- 
mal shock  is  propagating  is  characterized  by 

U1  = V1  = 0,0  ’ 

P1  = 1.01325x10s  [Pa]  , 

Pj  = 1.225570786  [kg/m3]  , 

(C3) 

e2  = 2.06689408xl05  [J/kg]  , 

where  e is  computed  from  the  perfect  gas  formula  e = 2.5p/p  for  y = 1. 
The  Ranxine-Hugoniot  relations  for  a normal  shock  wave  in  a perfect 
gas11  imply: 


(C4) 


Pj  (Y-l) 

P2(W)1 

2pl 

* J 

1 P2 

Y-l  P2 

> 

1/2 

> 


^H.W.  Liepman  and  A.  Roshko,  Elements  of  Gas  Dynamics,  Wiley  and  Sons , 
1967,  pp.  62-65.  ' 
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where  the  subscript  1 and  2 denote  the  value  of  the  variables  in  the 
quiescent  and  disturbed  states,  respectively,  cs  is  the  shock  front 
velocity  and  y = 1.4.  In  particular  for  a shock  of  strength  (P2/P1) 
five  and  a quiescent  state  characterized  by  equations  (C3) , we  have 

v 2 = 4.6190563x102  [m/s]  , 

u2  = 0.0  [m/s]  , 

p?  = 5.06625X105  [Pa]  , (C5) 

P2  = 3.4538813  [kg/m3] , 
e0  = 3.6670701xl05  [J/kg]  . 
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*♦♦♦  THE  FOLLOWING  PILOT  HYDRO-CODE  USING  A FI NI TE-ELEMENT 

METHOC  TO  SOLVE  THE  UNSTEADY  GAS  FLOW  EQUATIONS.  THE  FOLLOWING  LISTING 
CORRESPONCS  TO  THE  VERSION  OF  THE  CODE  WHICH  PRODUCED  FIGURES  7,8,9,10  OF 
THIS  8RL  REPORT 

CONTENTS  OF  COMMON  BLOCKS  BL OK  1 , BL0K2  AND  BL0K3  ARE 
XVI1,J)*ABSCISSA  OF  NODE  J 
XVI2,J»=0RDINATE  OF  NODE  J 

K»K0LC  OR  KNEW-T1ME  LEVEL  AT  WHICH  THE  ABOVE  VARIABLES  ARE  KNOWN 
OR  UNKNOWN 

VARIK, 1,J)*ABSCISSA  VELOCITY  COMPONENT  AT  NODE  J 
VARIK, 2, J»*OROINATE  VELOCITY  COMPONENT  AT  NODE  J 
VARIK, 3, J)«PRESSURE  AT  NODE  J 
VAR(K,A,J)«SPECIFIC  INTERNAL  ENERGY  AT  NODE  J 

VAPT IL,JI=PARTIAL  DERIVATIVE  OF  VAR  I KNE W ,L , J I WITH  RESPECT  TO  TIME 
AT  NOCE  J 

NR  I M , J I = NUMfi E R OF  THE  NEIGHBORING  NOOE  M TO  NODE  J 
DEFINITION  OF  VARIABLES 

NUMPER=TOT AL  NUMBER  OF  NODES  IN  THE  PROBLEM, .LE . 2000 

NUMONC  = TOT AL  NUMBER  OF  NEIGHBORS  FOR  A GIVEN  NODE 

NOEC=NUMPER  OF  DIFFERENTIAL  EQUATIONS  TO  BE  SOLVED,  MOST  OFTEN  A 

MU A = NUMBER  OF  UNKNOWN  PARAMETERS  All)  I =1 , 2 , . . . MU A , MOST  OFTEN  12 

ORGT I M=OR  IG IN AL  OR  INITIAL  TIME 

MXNTMS=MAX IMUM  NUMBER  OF  TIME  STEPS 

C0MM0N/BLCKl/XYI2,2CC0),NRI9,2C00)/BLOK2/VAR(2»A,200O)/BL0K3/VARTI 
1 A , 2000 ) 

FOR  SHOCK  PROBLEM,  THE  SHOCK  STRENGTH  IRATIO  OF  PRESSURE)  MUST  BE 
SFT  IN  BOTH  EXTSOL  ANC  GRAPHIT  ♦•MANUALLY*** 

CERTAIN  OTHER  PARAMETERS  MUST  BE  SET  MANUALLY  IN  EXTSOL 

CALL  STARTINUMBER,NUMOND,NOEQ,MUA,KCLD,KNEW,CRGTIM,MXNTMS,BEGPOS) 

INCEX*0 
OT IME  = 0RGT  IM 
T I MF*ORGT  IM 
19  INCEX= INCEX*  1 

CALI  T I ME STE PI  NUMBER, NUMOND,NCEQ,MUA,KOLO, KNEW, OTIME, TIME  .INDEX) 

CALL  CISPLAYINUMBER, KOLD, KNEW, TIME, INDEX, ORGTIM.BEGPOS) 

I F I INDEX  .LT.  MXNTMSI  GO  TC  19 
PRINT  52 

52  FORMAT! 1H.13X, 'PROGRAM  IS  F IN  I SHED- S TOP • ) 

STOP 

ENC 

SUBROUTINE  ST  ART ( NUMBER, NUMONC , NOEQ , MUA , KOLD , KNE W ,ORGTI M , MXNTMS , 
l BEGPOSI 

♦♦START-GENERATES  ENTRIES  FOR  ARRAYS  X Y,NR , VAR ( KNE W , A ,2 000 ) , V ART 
-ASSIGNS  NUMBERS  TO  ITS  ARGUMENT  LIST 
C0MMON/PLCKI/XYI2,2000),NR(5,200O)/BLOK2/VARI2,A,2COO)/BLOK3/VARTI 
1 a , 2000 ) 

C ♦♦  THE  FOLLOWING  ARE  GEOMETRIC  PARAMETERS  FOR  A NORMAL  SHOCK  PROBLEM 
C Y1=INITIAL  VALUE  OF  Y IN  DOMAIN 

Y 1 *2 .965 

C Y2 *F I N AL  VALUE  OF  Y IN  DOMAIN 

Y2*3 .050 

C C Y*S 1 2 E OF  MESH  IS  Y DIRECTION 

C Y *0 .00  5 

C N0XC*NUMPER  of  divisions  in  X-OIRECTION 
NOXC*A 

C eECPOS  IS  BEGINNING  POSITION  CF  SHOCK 

45 
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PRECEDING  page  BLANK- GOT  FILMED 


BEGPOS*  3.0025 

C XBEG I N* IN  IT  l AL  X POSITION 

X B EG  I N = 0 . 0 
ORGT I M=0 .0 
MXNTMS=14 

C*U  MESGEN(NUMBER»Y1,Y2»DY»NCXD»XBEGIN) 

PRINT  103 

103  FORMAT!  1 F , // » 20X  , * MATRICES  XYC2.2C00I  AND  NR(9,2000)  ARE*, //I 
PRINT  104 

104  FORMAT  < 1H , 5X , • NODE • . 13X , • CORRDINATES • ,25X , • NE I GHBOR I NG  NODESIl-8 
1)' ,24X,'N0CE  ') 

PRINT  105 

105  FORMAT  ( IF , 10 5X , • T YPE  • ) 

PRINT  101 

101  FORMAT  ( 1F.6X, • I* ,11X, 'X', 17X, 'Y', 13X, 'l',6X,'2' ,6X,'3* ,6X,'4' ,6X, 
1'5',6X,,6',6X,,7,,6X,,8',6X,,9') 

CO  51.  1=1, NUMBER 

51  PRINT  102,  I,XY!  1,  I ),XY( 2,  1 >,(NR( J,I ),J  = 1,9) 

102  FORMAT  (IF,  2X  , 1 5 , 3X , 1PE  1 5 . 8 , 3X  , IPE 1 5.  8 , IX  , 9(2X,15IJ 
NUM0NC=9 

N0EC=4 
KNEW=l 
K0LC=2 
MU  A = 1 2 

CO  10  N0NUM=1, NUMBER 
X*  X Y ( 1 , NONUM I 
Y*XY ( 2.N0NUM ) 

T=ORGTIM 

CALL  EXTSOLIX.Y.T , UV , VV , P , E , SPECV, IBADI 
IF(  IBAC.NE.O  IGOTO  20 
VAR (KNEW, 1 , NONUM ) *UV 
VA«(KNFW,2,NCNUM)=VV 
VAR(KNEW, 3 .NONUM ) =P 
VAR(KNEW,4,NGNUM)=E 
10  CONTINUE 

CALL  IN  I TNG ( NUMBER ,NUM ON D, NOEC, KNEW, MUA I 
CALL  OUTN04INUMBER,KNEW,ORGTIM) 

RETURN 

20  PRINT  21, NONUM 

21  FORMAT! 1H,/, • TROUBLE  WAS  INCURRED  IN  FINDING  EXACT  VALUE  AT 
10CE=  *,16,/,'  PROGRAM  IS  STOPPED  BY  PROGRAMMER') 

STOP 

RETURN 

ENC 

SUBROUTINE  T I MESTEP ( NUMBER , NUMONO.NOEQ ,MUA,KOLO .KNEW ,0T l ME ,T IME, 

1 INCEX) 

C COMPUTES  VALUES  OF  FLOW  VARIABLES  AT  INITIAL  TIME*DT  AT  ALL  NODES 
C KOLC  INCICATES  LEVEL  OF  KNOWN  INFC  IN  VAR 
C KNEW  INCICATES  LEVEL  OF  UNKNOWN  INFO  IN  VAR 

COMMON/BL0Kl/XY(2,200OI,NR(9,2CCO)/BLOK2/VAR(2,4,2O0O)/BL0K3/VART( 

14,2000) 

C 4*THE  FOUR  KEY  SUBROUTINES  OF  CALCI  ARF  TRANSL , SOL VER , 1 NI TZR , VALNEW 
EXTERNAL  TP ANSL , SOL VER , IN  I TZR , VALNE W 
C SWITCH  INDICATORS  FOR  CALCULATION  AT  NEW  TIME  LEVEL 
OT  IME=TIME 
I S AVE*KNEW 
KNEW=KOLC 
KOLC*  ISAVE 

C COMPUTE  TIME  INCREMENT 
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CALL  DELTAT I NUMBER, NUMOND, NOEQ, KOLO, OT ) 

C COMPUTE  NEW  TIME 
T I ME*OT I ME  ACT 
PRINT  4000,01  IME.CT  .TIME 

4000  FORMAT  ( IF , / , 7X , • SUB  TIMESTEP  OLD  T I ME* • . IPE 1 5 . 8 , / .22 X . • OELTAT* • 
1,1PE15.8,/,22X,  ‘TIME  TO  BE  USED  IN  CALCULATIONS*  IPE 15.81 
C LOOP  OVER  NOCES 

C ^BECAUSE  ALL  THE  VALUES  OF  THE  FLOW  VARIABLES  ARE  EQUAL  FOR  CONSTANT 
C V VALUES  ONLY  ONE  VALUE  AT  A GIVEN  V MUST  BE  CALCULATED 

NONUM=  3 
28  CONTINUE 

C IF  NR1NU MONO. NON UM)=l,US£  C ACC  I ( NON-BDY  NOOEI.ELSE  CALCBIBOY  NODE  I 
IF(NR(NUMCNC,NONUM).EQ.l»GOTO  25 
CALL  CALC0 ( NONUM. KNEW. TIME. INDEX) 

GOTO  20 

25  CALL  CALC  I ( NUMBER , NUMOND , NOEQ ,MU A ,KOLD, KNE W , T I ME ,0 T .NONUM , TRANSL , S 
iOLVFR. INITZR.VALNEWI 

C **USE  THE  VALUE  AT  THE  NODE I NONUM ) TO  COMPUTE  VAR  AND  VART  ENTRIES  FOR 
C A GIVEN  Y VALUE 

20  CHKRC=XY<2, NONUM) 

RACVEL=VAR(KNEW,2,N0NUM) 

N=N0NUM-3 

26  N = N ♦ i 

IF<  N .CT.  NUMBER)  RETURN 
RC=XYI2,N) 

I F ( R C .GE.  (CHKRDaO.OCOOOI ) ) GO  TO  27 
VAPIKNEW, 1, N ) =V AR ( KNEW, 1. NONUM) 

VARIKNEW.2.N  »=RADVEL 
VAR(KNEW,3,N)*VAR(KNEW,3«  NONUM I 
V API  KNEW, 4. N )=VAR( KNEW, 4, NONUM) 

VAPTI 1,N)*VART(  l, NONUM) 

VART(2,N)=VART(2,N0NUM) 

VAPT(3,N)*VART( 3, NONUM) 

VAPT(4,N)=VARTI4, NONUM) 

GO  TO  26 

27  N0NUM=N*2 
GO  TO  28 
ENC 

SUBROUTINE  C I SPL A Y ( NUMBER , K OLO ,KNE W , T I ME , I NOE  X ,ORGT I M .BEGPOS ) 

C ♦•CISPLAY-GENERATES  OUTPUT  AFTER  THE  FLOW  VARIABLES  ARE  COMPUTED  AT  A NEW 
C TIME 

C0MMON/BL0K2/VARI  2,4, 2000  )/BLCK  I/XY  ( 2 , 2000 ) , NR  I <3 ,2000  ) 

C **0L TN06  GIVES  LISTING  OF  VAR  AND  DENSITY  AT  NEW  TIME 
CALL  OUTN06INUMBER, TIME, KNEW) 

C ♦♦GRAPHIT  USES  THE  CALCOMP  PLOTTER  TO  GRAPH  PRESSURE , Y-VE LOC I TY , INTERNAL 
C ENERGY  PER  UNIT  MASS  AND  CENSITY 

CALL  GRAPFIT  I T IME .KNEW, NUMBER , INDEX , ORG T IM , BEGPOS ) 

RETURN 

ENC 

SUBROUT  INE  IN  I TNG! NUMBER .NUMOND , NOEQ .KNEW ,MUA ) 

C COMPUTES  INITIAL  VALUES  OF  THE  PARAMETERS  OF  FLOW  VARIABLES  AT  THE  INITIAL 
C TIME  LEVEL  ONLY 

COMMON/BL0Kl/XY(2,2000),NR(9,2C0O)/BLOK2/VAR(2,4,20O0)/6L0K3/VARTt 

14,20001 

Cl  MENS  ION  ALPHA(12),RH0INT(3,9),X(3),Y(3),GT(4,3I,G(4,3),D0GI4,2), 

1 A I I 12  ) , RHOT I 3 ),TRHOIN( 3, 3) 

C I M TNG  IS  WRITTEN  FOR  MUA*12.  CHANGES  MUST  BE  MAOE  IF  IM  TNG  IS  NOT  12.  INITNC 
C CHECKS  FOR  MUA-12. 
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IFIMUA.NC.12)G0T0  101 

THE  NUMBER  OF  NOCES  IN  OOM A IN  I NUMBER  I CANNOT  BE  GREATER  THAN  2000.  INITNG 
CHECKS  NUMBER 

IFINUNBER.GT.2COOIGOTO  201 
LOOP  OVFR  ALL  NOCES 

CO  20  NONUM*  1. NUMBER 
CO  25  J* 1 • NU * 

25  ALPHAI J ) =0.0 

COUNT  NUMBER  OF  NODES  PER  PRISM  — NDPPRM 
NC  PPRM  = 1 
NONMl-NUMONC-l 
CO  30  J* 1 .N0NM1 
IFINRIJ, NONUM). EO.OIGOTO  35 
30  NCPPRM=NCPPRM*1 

COMPUTE  VARIABLE  FROM  EQUATION  OF  STATE  AND  ITS  FIRST  DERIVATIVES  AT  PRISM 

NOOES 

35  CALL  INARES ( NUMOND, NOEQ, KNEW, NONUM.NDPPRM.RHO I NT ) 

TRANSFER  CENTER  NOCE  DATA  TO  WORKING  ARRAYS 
X ( 1 l=«XY(  l, NONUM  ) 

YI1MXYI2, NONUM) 

CO  40  1 1 1 .NOE  Q 
40  G( I. 1 )=VAR|KNEW,I,NONUM) 

DETERMINE  NUMBER  OF  TRIANGLES  PER  PRISM  NDPPRM-1  FOR  NON-BDY  AND  NDPPRM-2  FOR 
BTY  NOCE 

NUMTR I=NCPPRM-1 

IF  (NRINUMONC.NONUMI.NE.l  )NUMTRMN0PPRM-2 
LOOP  OVER  THE  TRIANGLES  OF  THE  PRISM 
RNMNUMTR I 
CO  45  I TR I * l . NUMTR I 
THIRC  VERTEX  OF  TRIANGLE  IS 
I V3=  I TR I ♦ 1 

EXCEPT  WHEN  CENTER  NOCE  IS  NOT  ON  BOUNDARY  AND  LAST  NGHB.  NODE  OF  PRISM  IS 
CCNSICEREC 

IFIITRI.EC.NCPPRM-l)IV3=i 

TRANSFER  DATA  OF  NEIGHBORING  NOOES  TO  WORKING  ARRAYS 
JSUP2*NR ( ITR I .NONUM ) 

JSUR 3 = NR ( IV3, NONUM  I 
X ( 2 > «XY ( 1 . JSUB2  ) 

X ( 3 ) =XY ( 1 . JSUB3  ) 

Y12)=XY(2»JSUB2I 
YI3MXYI  2.JSUB3) 

ABSOLUTE  VALUE  OF  CET  IS  THE  AREA  OF  THE  TRIANGLE/2 

CET  = XI1MYI3)-XI3)*Y(1MX(2)MV(1)-XC1)*YI2MXI3MV|2)-XI2)*Y(3I 
IFICET.EC.O. ) GOTO  19 

DOGI I»J)  = PARTIAL  DERIVATIVE  OF  VAR  I KNEW . I .NONUM I WITH  RESPECT  TO  XIJ=1)  OR 
Y I J=  2 ) ON  THE  TRIANGLE  WITH  CENTER  AT  NONUM 
CO  50  M 1 »NOEQ 
C< 1,2 )=VAR( KNEW, I , JSUB2 ) 

G ( I,3MVAR(KNEW,I,JSUB3) 

COG ( 1,1  Ml  (G( 1,1 )-G  I 1,2) )*YI3I-IGII»1)-G(I»3) ) * Y 1 2 ) ♦ IG  1 1 ,2 l-GI I ,3 
1 I MYIi) l/OET 

50  COGI I » 2 M I - 1 G I I , 1 )-G I 1,2) )*X(3M(G(I,1)-G(I,3) )*XI2)-(Gll»2l -G (1,3 
1 I ) *X ( 1 ) l/CET 

TRANFER  INFO  FROM  RHOINT  TO  TRHOIN  AND  COMPUTE  APPROXIMATE  DERIVATIVES  OF  U,V 
AND  INT  ENG  WITH  RESPECT  TO  TIME  BY  SOLVING  FOR  TIME  DERIVATIVES  IN  THE 
GOVERNING  EQUATIONS 
CO  55  Ml, 3 

TRHOINI  I,  IMRH0INT1  1,11 
TRHOINI  I.2MRH0INTI  I, ITR I «1 I 
TRHOINI  1,3  MRHO  INTI  I , I V 3 ♦ 1 ) 


A 
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GT (l, ll*-!G!l, I KCOGI 1,1I'G(2,I  >*D0G<1, 2>*DOGI3,l)/TRH01N(l,l  » I 
GT12.I >»-!G!  l,  I 1 *COG  ( 2 , 1 I *G ( 2 » I l*DOG  I 2, 2 I »DOG  1 3 , 2 1 /TRHO  IN  1 1 , I I I 
55  GT 14, I l*-!G! 1, I )*DOG( A, 1 DGI  2,1 )*DOG 1 4, 2 1 4G I 3 , 1 I * (DOG (1,1 ) ♦DOG (2 >2 
1)I/TRH01N!1,I>> 

C CONFUTE  INITIAL  VALUES  OF  PARAMETERS  FOR  U.V,  INT.  ENG. IAfl-6. 10-121  I 
C SOLVE  FOR  A< 1I,A(2),A(3I  BV  EQUATING  A ! 1 1 • X 1 1 I ♦ A 1 2 > * Y ( I 1 ♦ A I 3 1 *G  T 11  , 1 1 , I = 1 , 2 , 3 
CO  60  I*i,NOEQ 
IFI I.EQ.3IGOTO  60 
J*  I * 3 
JM1=J-1 
JM2= J-2 

Al(J)=<GT(I,l )*!X! 3>*Y!2I-X(2)*Y! 3 1 1 +G  T II  , 2 1 * ( X ! 1 1 * Y 1 3 1 - Y 1 1 1 *X ! 3 1 1 
1 »GT ( I,3)*(X(2I*Y( 1 1-X I 1 1 <Y I 2 I ll/DET 

A I ( JMi  I * ( -GT  ( 1,1 >*<X<  3I-XI21  l«GTII.2l*(Xm-Xlll  1-GT 1 1 , 31 ♦ I X 12 l-X ! 

III ) I/DET 

AHJM2»*(GT!l,ll*<Y!3)-Y!2>)-GT!I,21*IYI3l-YIll  1 +GT II ,3 1* I VI2 1-V 1 1 
1 ) I l/OET 
60  CONTINUE 

C CONFUTE  DERIVATIVES  OF  DENSITY  ANC  PRESSURE  WITH  RESPECT  TO  TIME  ANO 
C PARAMETERS  A ( 7-9  > FOR  THE  PRESSURE 
CO  65  1*1,3 

C KHOT  IS  THE  TIME  DERIVATIVE  OF  DENSITY 

RHOT ( I )*-(G( 1, I )•( TRHOINI 2, I 1 <DOG! 3, 1 1 ♦TRHO I N I 3 , I I *DOGI 4, 1 1 I 

1 +G(2, I I<<TRHOIN< 2, I )*DOG!3,2DTRHOIN!3,ll*OOG!4,21) 

2 ♦TRHO IN < 1, I 1* ( OCG! 1,1) ♦DOG (2, 21  1 1 

C GT ( 3 , I I IS  THE  TIME  DERIVATIVE  OF  PRESSURE  AT  VERTICES 

6 5 GT  < 3 , I ) = ( RHOT I I l-TRHO IN < 3, I I * < A I < 10  I *X < I) ♦A  I < 1 1 I * Y ( I ) ♦ A I < 12  I I) /TRH 
10  I M 2,  I I 

AI!7I=!GT(3,1)*!YI3I-Y!2) I-GT < 3 ,2 ) * < Y < 3 )-Y < 1 11 ♦GT 1 3 ,3 1 * I Y 1 2 l-YI 1 ) I 
1 I/CFT 

A I (8 ) = <-GT<  3, 1 1*1X1  3)-X<  2)  I +GT I 3, 2 ) * < X I 3 » -X < 1 1 1 -GT < 3 , 3 I * < X < 2 1 -X (II 
1 I I/DET 

Al ! 9 1 = 1 GT (3,1 )*(XI31*Y!21-X(21*Y!31 •♦GTI3»21*!XI1I*VI3I-Y!1I*X13I I 
l ♦GT!3,31*(X!2I*YI ll-XI 1 )*Y! 2 1 1 1 /DET 

C SUM  THE  tlll'S  OVER  THE  TRIANGLES  IN  THE  GIVEN  PRISM 
CO  70  1*1, MUA 

70  ALPHA!  I l*ALPFA(  I MAH  l I 
45  CONTINUE 

C THE  AVERAGE  OF  THE  AID'S  IS  THE  INITIAL  VALUE  OF  THE  A(I1'S  FOR  THE  PRISM 
CO  75  1*1, MUA 
75  Aid  1 * AL  PH  A ( I I/RNT 

C STORE  IN  MATRIX  VART  THE  INITIAL  VALUES  OF  THE  PARAMETERS  AT  THE  CENTER  NODE 
C IN  THE  TRANSLATED  COORDINATE  SYSTEM.  WE  TAKE  A I 1 , 2 , 4 , 5 , 7, 8, 10  AND  111=0 
C ANC  A(  I*3»  = VART(  I , NONUM I 
CO  AO  I * 1 ,NOEQ 

PO  VART<  I , NONUM  I *A  II  1 * 3-2 1 *XY ( 1 , NONUM 1 ♦ A I ( I* 3- 1 1* X Y ( 2 .NONUM I ♦ A I d *3 1 

20  CONTINUE 
RETURN 

101  PRINT  102, MUA 
STOP 

201  PRINT  202, NUMBER 
STOP 

IV  PRINT  5, NONUM 
STOP 

3  FORMAT! IH,//, • CET*0.0  IN  SUB.  INITNG',/,  • CNE  CF  THE  TRIANGLES 
1 ASSOC  I AT  EC  WITH  NODE', 15,'  HAS  ZERO  AREA',/,'  CHECK  THE  ORDERING  A 
INC  NUMBERING  OF  NGHB.  NODES',/,'  STOP'I 

102  FORMAT  ! IF , • MUA*', 15,'  BUT  SUBROUTINE  INITNG  IS  WRITTEN  FOR  MUA* 

112, STOP* 1 

202  FORMAT  I1H,  • NUMBER* • , 1 7, • BUT  INITNG  IS  WRITTEN  FOR  NUMBER  NOT 
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1GREATER  THAN  2CQ0.  STOP*  I 
ENC 


SUBROUTINE  MESGFN ( NONUMF ,R 1 ,R2, OY.NOTO, THE TA8 I 
THIS  SET  OF  SUBROUTINES  GENERATE  A FINITE  ELEMENT  MESH  IN  CARTESIAN 
COORCINATESCX.V)  THEY  WERE  ADAPTED  FROM  THE  CLYDRICAL  COOR.  PROGRAM 
R IS  THE  Y COORDINATE  AND  THETA  THE  X COORDINATE 

•♦FINITE  ELEMENT  GENERATING  PACKAGE  MAY  BE  CONFUSING  IN  ITS  LOGIC 
BECAUSE  IT  HAS  BEEN  ADAPTED  FROM  POLAR  TO  RECTANGULAR 
GEOMETRY 

COMMON/BL0K1/XYI2.200OI.NRI 9.20C0I 
C I PENSION  RAC  IK  200  l 
C NOTD  IS  NUMBER  OF  DIVISIONS  IN  THE  X DIRECTION 
NOTCP1=NOTC*1 

C THETAB  IS  INITIAL  X POSITION 
C CY  IS  INCREMENT  IN  Y-CIRECTION 
C R 1 IS  INITIAL  VALUE  OF  Y 
C R2  IS  FINAL  VALUE  OF  Y 

C OR  IS  INCREMENT  IN  X-CIRECTION  TO  INSURE  EQUILATERAL  TRIANGLES 
CR*l.i54700538*CY 
THETAE  = THETAB'»NOTD*DR 

C NOCPOS  INDICATES  TYPE  OF  POSITION  THE  FIRST  NODE  WILL  HAVE-SEE  AUTGEN 
N0CP0S=1 
NONUMB*NOT0+2 
RACIim-RI 
1*1 

23  I>IM 

RACIIII)*RADII( I-1I+0Y 

i f i (RAoiim*o.ocoon  .lt.  R2i  go  to  23 

PRINT  24  » (RAC  III IKL  I » IKL“1»  I I 

24  FORMAT  <1F,10X,'TFE  Y SUBDIVISIONS  FROM  SUBPROGRAM  MESGEN  FOR  X-Y 
1 FIN.  EL.  MESH  ARE • . I / , 1 3X , 1PE 1 5 . 8 II 

C AUTGEN  GENERATES  GRID  FOR  Y.GT.Y1-R1 

CALL  AUTGEN  ( R AD  I I , I , NODPCS.NCNUMB, NOTD, THET AB . THET AE .NONUMF I 
C GENERATE  GRIC  FOR  Y*Y1*R1 
CO  20  I * 1 .NOTOP 1 
NR  1 9 » 11*7 
CO  20  J*  5»  8 

20  NR  I J.  1 1*0 

CO  21  I * 2 »NOTCP 1 
NR (It!  1*1-1 
NR ( 2 * I I * I ♦NOTD*  1 
NRI 3»  I I*I4N0TD«2 

21  NR(4.II>I*1 
NRIl»l»*N0TC+2 
NRI2.  ll-NCTCO 
NRI3, 11*2 

N R ( 4 » 1 1 *0 
NRI4. 161*0 

CTFETA=( THE TAE-THE TAB! /NOTD 
CO  22  I * 1 , NOTCP 1 
R I * I -1 

THETA*THETAE+OTHETA*RI 
X Y ( 1 • II  * THET  A 

22  XV(2»II*R1 
RETURN 
ENC 


SO 


SUBROUTINE  AUTGENIRADI I , N .NODPOS.NONUMB, NOTD , THE  TAB , THE TAE , NONUMF I 

C THE  AUTOMATIC  GENERATION  OF  THE  FINITE  ELEMENT  GRID  STARTS  WITH  R!« 

C CELTAR1.RACIII2) 

C NOT C IS  NUMBER  OF  THETA  DIVISIONS,  NONUMF  *TOT AL  NUMBER  OF  NODES  INONUMB* 

C THOSE  GENERATED 

COMMON/BLOK1/XY(2,2COO),NR(9,20CO> 

C I MENS  ION  RAC  1 1 1 200  I 
PI*3. 1A159  26  5 35897932A 

C NONUMB=THE  BEGINNING  NODE  NUMBER  AT  RADIIC2). 

NONUM»NONUMe-I 
NDP1*N0T0'*  1$  NCP2*N0TC+2 
CTHFT  A*  < THETAE-THETABI/NOTD 

C N0CP0S=0  IF  NOCE  NONUMB  IS  A VERTEX  OF  A COMPLETE  ECLILATERAL  TRIANGLE, 

C IF  NOT  NODPOS*  1 
KNOW*NODPOS 
KK.KNOW* I 

KBEF0R*XMCDF(KK,2) 

CO  20  I»2,N 
IFIKNOW.NE.l  I GOTO  60 

C VERTEX  OF  AN  INCOMPLETE  EQUILATERAL  TRIANGLE  IS  ON  THE  RIGHT  BOUNDARY 
IF ( I . EQ.N (GOTO  30 

C RI.LT.RACI I .LT.R2  AND  THETA.THETAB 
NONUM-NONLMM 
XYU,N0NUM)=TFETA8 
XYI2, NONUM). RACIII I > 

NR19, NONUM )=2 
NR(1, NONUM). N0NUM-NDP1 
NR  12  .NONUM  » = KONUM«i 
NR( 3.N0NUM ).N0NUM«N0P2 
CO  36  J.  A , 8 
36  NR ( J , NONUM  ) *0 

THETA*THETAB*0.5*CTHETA 

C R1.LT.RACU.LT.R2  ANC  THETA*THET AB 

C Rl.LT.RACIl.LT.P2  AND  THETAB .L T . THETA .L T. THE TAE 
AO  N0NUM*N0NUM»1 

XYI1, NONUM). THETA 
X Y ( 2 , NONUM ) .R  AC  1 1 I I I 
NR(9, NONUM). 1 
NRU, NONUM). NONUM-1 
NR(2,  NONUM)  .N0NUM-*NCP2-l 
NR (3,NONUM)*NONUM»NCP2 
NRIA, NONUM). N0NUM«1 
NRI5,N0NUM).N0NUM-NCP1 
NR(6, NONUM). NONUM-NOPI-l 
NRIT, NONUM). C 
NR(8, NONUM). 0 
THETA. THETA»CTHETA 
IFITHETA.LE.THETAE )GOTO  AC 

C Rl.LT. RADII. LT.R2  AND  THET A*THET AE 
NONUM*NONUM* 1 
XYI1, NONUM). THETAE 
XYI2, NONUM). RADIIIII 
NR ( 9, NONUM  » ■ 3 
NRI 1, NONUM). NONUM -NDP2 
NRI2, NONUM). NONUM-1 
NR(3,N0NUM)«N0NUM*NCP1 
CO  39  J. A , 8 
39  NRI J, NONUM). 0 
GOTO  50 

C VERTEX  OF  A COMPLETE  EQUILATERAL  TRIANGLE  IS  ON  THE  RIGHT  BOUNDARY 

SI 


60  IF  ( I . EQ.N  IG0T0  35 
N0NUP-N0NUP41 

C Rl.LT.RAClI.LT.R2  AND  THETA-THETAB 
XYU,NONUP)»THETAB 
XY<2,N0NUP>»RACIII  I ) 

NR  1 9 .NONUP  > *2 

NR  1 1 ,N0NUP)»N0NUM-NCP2 

NRC2,NONUPI»NONUM-NOP24l 

NR  I 3 .NONUP  ) »NONUM*  l 

NRIA.NONUPI-NONUMaNDPI*! 

NR  C 5 .NONUP l*N0NUM4NDPl 
NR ( 6 • NONUP 1*0 
NR ( 7«N0NUP l»0 
NR  ( 8 . NONUP 1*0 
THETA^THETAB+DTHETA 

C Rl.LT. RADII. LT.R2  AND  THETAB .LT . THETA .LT. THE TAE 

A A N0NUM»N0NUP*1 

XVI  1,N0NUP)*TFETA 
XY 12,N0NUP)*RAD1I I I > 

NR ( 9. NONUP )*l 

NR( 1 . NONUP )*NONUM- 1 

NR  ( 2 .NONUP I *NONUP*NCP 1 

NRI3,N0NUP>«N0NUM+NCPIU 

NR(A,NONUPI»NONUM4l 

NRI5.N0NUP  »*N0NUM-N0P24l 

NR(6,N0NUP)*N0NUM-NCP2 

NR ( 7 , NONUP 1*0 

NR  I 8 .NONUP ) *0 

THETA*THETA+CTHETA 

IFITHETA.LT.(THETAE-0THETA/2. » IGOTO  AA 

C Rl.LT. RACI1.LT.R2  AND  THETA  = THETAE 
N0NUP*N0NUP4l 
XYU.NONUP»  = TFETAE 
XYI2.N0NUP >*RAD 1 1 ( I ) 

NR  1 9. NONUP ) * 3 
NR(l.NONUP)*N0NUM-NCPl 
NR ( 2 » NONUP ) *NONUP-NOP 1-1 
NR ( 3 « NONUP I * NONUH- 1 
NR(A.N0NUP»«N0NUM4N0P2-1 
NR  1 5, NONUP I«N0NUM*N0P2 
NR  ( 6 . NONUP > *0 
NR  1 7, NONUP >«0 
NR  I 8 * NONUP ) *0 

50  I $ AVE^KNON 
KNOH'KBEFOR 
K BEFOR= ISAVE 
GOTO  20 

30  N0NUP-N0NUP41 

C RACII-R2  ANO  THETA-THETAB 
XYIl.NONUPI-THETAB 
XY(2,N0NUP»-RACII(I) 

NR (9 .NONUP ) *5 
NR  1 1 » NONUP ) -NONUH 4 1 
NR  1 2 . NONUP I -NONUM-NDP 1 
CO  51  J«3,8 

51  NR  I J » NONUP I »0 
THETA»THETAB*0.5*CTHETA 

62  NONUH-NONUH* 1 

C P AC  1 1 *R2  and  thetab.lt. theta. lt.thetae 

XYI1.N0NUPI-THETA 


52 


XV ( 2 • NONUN I *R AC  1 1 1 l I 
NR(9, NONUNI-4 
NR  U, NONUN  »*NONUM*l 
NRI2, NONUNI -NONUM-NDP 1 
NR ( 3 ,NONUN > *NONUH-NDPl-l 
NR (4 • NONUN l-NONUM- 1 
CO  52  J*5,8 

52  NR  ( J , NONUN  1*0 
THETA«THETA*CTHETA 
IFITHETA.LE.THETAE IGOTO  62 
NONUH-NONUMM 

C RACII-R2  AND  THETA-THETAE 
XV ( It  NONUN l-THETAE 
XV(2*N0NUN)=RACIt ( I ) 

NR 19, NONUN  1*6 

NR I 1 , NONUN  > -NONUM-NOP 2 

NR 1 2, NONUN l-NONUM-! 

CO  53  J=J,8 

53  NR ( J. NONUN ) *C 
GOTO  20 

35  NO NUN* NO NUN* 1 
C R AC  1 1 = R2  AND  THETA-THETAB 
XV(l,NONUNI*THETAe 
X V I 2 * NONUN ) *R AC  1 1 ( I ) 

NR I 9 »NONUN  ) * 5 
NR  1 1 f NONUN ) *NONUM«  1 
NR(2,NONUN)*NONUM-NCP2»l 
NR (3»NONUN)=NQNUM-NCP2 
CO  54  J = 4,8 

54  NR ( J ( NONUN  ) *0 

THET A-THETAB+CTHETA 
49  NONUM-NONUM*  l 

C RACI l =R  2 ANC  THETAR.LT. THETA.LT. TNETAE 
X V ( 1 , NO NUN  J * THET  A 
XV(2, NONUN I-RACIII I I 
NR  I 9 , NONUN 1 = 4 
NR(1,N0NUN)»N0NUN*1 
NR (2, NONUN  I -NONUM-NOP 2» 1 
NR ( 3,NONUN)*NONUM-NCP2 
NR  ( 4 , NO NUN ) * NON UN- 1 
CO  55  J-  5, 8 

55  NR(  I , NONUN ) *0 
THETA-THETA*CTHETA 

IFITHETA.LT.ITHETAE-OTHETA/2. I I GOTO  49 
NONUM-NONUM*  1 

C RACI I *R2  ANC  THETA-THETAE 
XVII, NONUN  l*THETA£ 

X V ( 2 , NONUN )*RACIII I > 

NR (9, NONUN | *6 

NR  I 1 » NO NUN  I *NONUM-NCP  l 

NRI2,N0NUN)*N0NUM-NCPl-l 

NR I3,N0NUN)*N0NUM-1 

CO  56  J * 4 , 8 

56  NRIJ, NONUN  1*0 
20  CONTINUE 

NONUNF-NONUN 

RETURN 

ENC 

SUBROUTINE  IN ARES (NUNONO ,NOEQ ,KNEN ,NCNUN .NOPPRN  ,RH01 NT ) 

S3 


C COMPUTES  VARIABLE  FROM  EQUATION  OF  STATE  AND  ITS  FIRST  DERIVATIVES  AT  PRISM 

C NODES 

COMMON/BLOK17XY(2,2000),NR(9,2000I/BLOK2/VAR  (2 ,*,20001 
C I MENS  ION  RHOINTI 3.NUM0ND > ,R I 3 I 

CALL  INEQSTIVARIKNEW,3,N0NUMI , VAR  I KNEW, A.NONUM ) ,R) 

CO  20  I*  l » 3 
20  RHO INTII,1)-Rt  I) 

CO  25  K-2.NCPPRM 
J SUP -NR  I K-l , NONUM  I 

CALL  INEQSTIVAR|KNEW,3,JSUB),VARIKNEW,A,JSUB>,R) 

CO  25  1-1,3 
25  RHOINTI I ,K )«R( I I 

C IF  NUMBER  OF  TRIANGLES  IS  LESS  THAN  MAX  NUMBER  ALLOWEO  IN  A PRISM,  SET 
C REMAINCFR  OF  RHCINT  TO  ZERO. 

I F I NOPPRM .EQ.NUMOND  I RETURN 
NFPP1-NCPPRMM 
CO  30  J-NPPP  1 , NUMCNO 
CO  30  1-1,3 
30  RHOINTI I, JI-0.0 
RETURN 
ENC 

SUBROUTINE  INECST I P, E , R ) 

C COMPUTES  THE  CENSITY  AND  ITS  CERIVATIVE  WITH  RESPECT  TO  PRESSURE  AND  INT. 

C ENERGY  FOR  A NOBEL-ABEL  GAS. 

C GAMMA  AND  ETA  ARE  TWO  PARAMETERS  IN  THE  NOBEL-ABEL  EQUATION  OF  STATE 
C TWO  PARAMETERS  ARE  GAMMA  ANO  ETA  IUNITS  FOR  ETA  ARE  M**3/KG» 

DIMENSION  R I 3 ) 

CATA  GAMMA, ETA/1. AO, 0.0/ 

G-GAMMA-1.0 

C-G*E«ETA*P 

CC-CAC 

C Klll-CENSITY 

r i n-p/c 

C RI2I-FIPST  CERIVATIVE  OF  Rill  WRT  PRESSURE 
RI2I«G*E/CC 

C RI31-FIPST  DERIVATIVE  OF  Rill  WRT  I NTFRNAL  FNERGY 
R I 3 I --GAP/CC 
RETURN 
ENC 

SUBROUTINE  OUTNOA I NUMBER , KNEW, CRGT I M » 

C *♦  OUTNOA-L I ST  S THE  VALUE  OF  MATRICES  VAR  AT  KNEW  LEVEL  AND  VART 
COMMON/BLOK2/VARI 2 , A, 2000  ) /BLCK 3/VARTI A, 2000  I 
PRINT  399 ,ORGT 1M 

399  FORMAT  I IF , / , IX , • THE  VALUES  OF  U,V,P,E  AND  RHO  AT  THE  INITIAL  •, 

1 ‘TIME  I ' , 1PEI6.8, * I ARE • , / I 
PRINT  397 

397  FORMAT! IH.5X, •NODE', 10X, 'U • , lex , • V • , 18X , 'P • , 1BX , • E • , 1 7X , • RHO* ) 

CO  II  1*1, NUMBER 

RHO- VAR  I KNEW, 3, I) / 1 0 . A*VAR I KNEW , A, I I ) 

11  PRINT  A02,I, IVARIKNEW.J, I ),J-1,A)»RHC 
A02  FORMAT! IH.3X, 1 5 , 5 1 AX , 1PE 1 5. 8 I) 

PRINT  398.0RGTIM 

PRINT  AOS, I I, I VART I J,  I ) , J - 1 , A ) , I - 1 , NUMBER  I 
39B  F0RMATI1H,/'  INITIAL  VALUES  CF  D/DT  OF  U,V,P,E  ARE  AT  TIME-', 

11PE15.B,//,5X,'N0CE',BX, • CU/DT • , 1 AX , 'DV/DT • , l AX , 'OP/DT' ,IAX,'DE/PT 

2'  > 

AC  a F0RMATI1H,I2X,I6,AIAX,1PE15.8)  M 
RETURN 
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ENC 

SUBROUTINE  EX TSOL < X , V, T, UV, VV,P ,E , SPEC V# I BAD » 

**  EXTSOL-FOR  A GIVEN  X,Y,T,EXTSOL  COHPUTES  THE  X AND  V VELOCITY  COMPONENTS, 
THE  PRESSURE, INTERNAL  ENERGY  BEHIND  A NORMAL  SHOCK  PROROGATING  INTO 
A CU I E SCENT  FIELD.  PERFECT  GAS  IS  ASSUMED  ON  EITHER  SIDE 
♦ ♦ VEL1,PI,R1,TI  ARE  Y-VELOC I TV .PRESSURE .OENSI TY , TE MPERATURE  OF  QUIESCENT 
STATE 

CAT  A VELI.P1.R1.T l/O.O, 1 . C 1 325E 5. I . 2255 T0786 .268 . 15/ 

C M ORGSHP,ORGSTM,SHKSTH  ARE  INITIAL  POS I TI ON. T I ME , STRENGTH  OF  SHOCK 
C M G IS  RATIO  OF  SPECIFIC  HEATS 

CAT  A ORGSHP.ORGSTM.G/3. 0025,0. 0.1. 4/ 

CATA  SHKSTH/  5.0/ 

C ♦♦  SPECV  ANO  IBAD  ARE  UNUSED  PARAMETERS 
SPECV-1.0  t IBAC-0 
C ♦♦  A1*S0UNC  SPEED 

AI=»SQRTIG*P1/R1) 

C0M=IG*l.0)/(G-1.0) 

C **  CS=SPEEC  CF  SHOCK  FRONT 

CS»A1*SQRT(  (G-l.O!/(  2.0*G)*(GU.0I*SHKSTH/I2.0*G)  ! 

POS  = V 

SHKPOS*ORGSHP+CS*( T-ORGSTM! 

IF  (POS  .GT.  SHKPOSI  GO  TO  15 

UP»( Al/G!*<SFKSTH-l.O)*SQRT  I ( ( 2 . 0*G I / ( G* 1 . 0 I ) / ( SHKSTH* 1. O/COM I I 

UV*0.0 

VV»UP 

P=SHKSTH*P1 

R = R1*(1.0«C0M*SHKSTH  )/(SHKSTH*COM) 

E *2 . 5*P/R 
RETURN 
15  UV*0.0 
VV*0.0 
P*P1 

E*2 . 5*P/R 1 

RETURN 

ENC 

SUBROUTINE  CELT  AT ( NUMBER ,NUMOND,NOEQ,K OLD, DT I 
C COMPUTES  THE  TIME  INCREMENT  FOR  THE  ENTIRE  DOMAIN  BY  USING  BI-CHAR.  ANO  VALUES 
C OF  THE  FLOW  VARIABLES  AT  TIME«OLD  TIME 

C0MM0N/PL0K1/XY(2,2CC0I,NR(S,2C00)/BL0K2/VAR(2,4,2000) 

INTEGER  START, SKIP, START2 

C START  IS  STARTING  NODE  SKIP  IS  NUMBER  OF  NODES  TO  BE  SKIPPED 
DATA  START.SK  IP/1,  1/ 

C FINC  MINIMUM  LENGTH(RADII)  OF  THE  FIRST  PRISM  ANO  SOUND  SPEED  AT  START  NODE 
J$UP*NR( 1, START) 

X R « ( XY( 1 , JSUB )-XY( l, ST ART) 1**2 

YR=(XY(2,JSUB)-XY(2, START!  )**2 

RM  IN  = XR*YR 

NONMl»NUMONC-i 

DO  20  I *2 , N0NM1 

JSUP*NR ( I, START) 

C IF  NR( I, START  )«0,  THEN  THERE  IS  NO  MORE  NGHB.  NODES  TO  THE  START  NODE, ELSE 
C THERE  ARE 

IF ( JSUB  .EQ.  0)  GO  TO  25 
XR * ( X Y ( 1 , JSUB  l-XYI 1, START  ) 1**2 
YR*(XV(2,JSUBI-XY(2,START)  )**2 
R * XR ♦ YR 

20  RMIN-MIN1FIRMIN.R  ) 

25  RMtN=SQRTF(RMIN) 
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CALL  SNDSPOIVARIKCLO, 3, STAR T » , VAR < KOLD , A , START ) ,A> 

C I NCLUCE  POSSIBILITY  FOR  UNKNOWN  (SWIRL  VELOCITY  IN  PCLAR  COORD) 

SWRLSP-0.0 

IFINOEQ.EC.5 ISWRLSP* VAR  I KOLD, 5, START) 

C FINC  THE  FIRST  VALUE  OE  DELTA  T 

TMIN-RMIN/ISCRTFI VARIKOLO,1,START)**2»VAR(KOLO,2,START)**2 
1 +SWRLSP*SWRLSP  )-»A) 

C INCREMENT  START  BY  SKIP  AND  PROCEED  THROUGH  THE  FINITE  ELEMENT  NET 
START2=START«SKIP 
CO  30  NONUM* ST  ART 2, NUMBER*  SKIP 
JSUP-NRI 1, NONUM) 

XR-IXYU,  JSUB»-XY(  l.NONUM) )**2 

YR-IXYI2, JSUB  )-XY I 2 * NONUM ) 1**2 

RM  IN*XR*YR 

NONMl*NUMONC-l 

DO  35  I*2,N0NMl 

JSUP  = NR(  1 1 NONUM  ) 

C IF  NR  I I • NONUM  )*0 • THEN  THERE  IS  NO  MORE  NGHB.  NOOES  TO  THE  START  NODE, ELSE 
C THERE  ARE 

IFI JSUB.EC.O (GOTO  AO 

XR-IXYI 1 « JSUB )-XY( 1, NONUM  ) )**2 

YR*IXV<2, JSUP  )-XY I 2, NONUM  ) )**2 

R*XR»YR 

35  RMIN*MINIF(RMIN,R) 

AO  RMIN*SQRTF(RMIN) 

CALL  SNDSPDI V AR I KOLD, 3 , NONUM ) , VAR  I KOLD, 4, NONUM ) , A) 

C INCLUDE  POSSIBILITY  FOR  UNKNOWN  (SWIRL  VELOCITY  IN  PCLAR  COORD  I 
SNRLSP=0 .0 

IFINOEQ.EC.5  ISWRLSP* VAR (KCLD, 5, NONUM) 

TM-RMIN/ISQRTFI VAR (KOLD, 1 .NONUM  )**2*VAR I KOLO , 2 , NONUM ) **2  + 

1SWPLSP*SWRLSP l+A) 

C FINAL  VALUE  OF  CELTA  T IS  THIN 
30  TMIN=MIN1F(TMIN,TMI 
CT*TM  IN 
RETURN 
ENC 

SUBROUTINE  CALC  1 1 NUMBER, NUMOND , NOE Q , MU A , KOLD .KNEW, TIME, DT, 

L NONUM, TRANSL, SOLVER, INI T2R, VALNEW) 

C COMPUTES  FOR  A NON-BOUNDARY  NODE  NEW  INITIAL  VALUES  FOR  THE  PARAMETERS  AND 
C VALUES  OF  THE  FLOW  VARIABLES  AT  THE  OLD  TIME  ♦ OT  (KNEW)  LEVEL 
EXTERNAL  TRANSL , SOLVER , IN  I TZR  , VALNEW, FU.FUSPEC 
C0MM0N/BL0K1/XY)2«2000),NR(9,2CC0)/BL0K2/VAR(2,4,2G00) 

1 /BLOK3/VART(4,2COO) 

C THE  VECTOR  OF  UNKNOWN  PARAMETERS  A IS  OFTEN  REDEFINEC  AS  B 

C THE  VECTORS  X AND  Y ARE  XY 1 1 , J ) , XY I 2 , J ) FOR  J CORRESPONDING  TO  NEIGHBORING  AND 
C CENTER  NOCE  OF  NODE  NUMBER  I NONUM ) . 

C MATRIX  VALS  CONTAIN  VAR  OATA  FOR  NEIGHBORING  NODES 
Cl  MENS  ION  B( 12),X(9),Y(9),VALS(A,9} 

C NCPPRM*NUMBER  OF  NODES  PER  PRISM  CENTERED  AT  NONUM  I.LE.9) 

CALL  TRANSL(NUMONC,NOEQ,MUA,KOLC,NONUM,NDPPRM,X,Y,VALS,B) 

CALL  SOLVER  I NUMOND.NOEQ, MUA , T IMF ,0T .NDPPRM, NONUM  ,X , Y , VALS, B,EU,FUS 
IPEC) 

CALL  INITZRINOEQ, MUA, NONUM, B) 

CALL  VAL NEW  I NUMOND, NOEQ, MUA, KNEW, DT, NONUM, VALS,B> 

RETURN 

ENC 

SUBROUTINE  CALCB I NONUM, KNEW, TIME, INDEX) 

COMMON/BL0Kl/XY(2,20C0),NR(9,20OO)/BLOK2/VAR(2,A ,2000 )/BL0K3/VART( 
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14,20001 
X»XY(  l, NONUM) 

Y = XY  I 2 .NONUM  I 
T*T I ME 

CALL  EXTSOL(X,Y,T,UV,VV,P,E,SPECV, IBAO> 

IF<  IBAC.NE.OJGOTO  20 
VARIKNEW, 1,N0NUM)=UV 
VAR  I KNEW ,2»N0NUM)*VV 
V AR ( KNEW , 3 , NO MUM ) *P 
VAR  I KNEW , 4 , NONUM I =E 
RETURN 

20  PRIM  2 1 , NONUM 

21  FORMAT! IE,/, * TROUBLE  WAS  INCURRED  IN  FINDING  BOUNDARY  VALUE  AT  N 
10CE=  • , 16,/,  • PROGRAM  IS  STOPPED  BY  PROGRAMMER*) 

STOP 

ENC 

SUBROUTINE  TR ANSL ( NUMONO , NOEQ ,MUA,KOLO, NONUM ,NOPPRM , X , Y , V ALS ,6  ) 

COMMON/ BLOK 1/XY I 2,20CO),NR(<9,2CCOI/BLOK2/VARC2,4,2COO)/BLOK3/VARTI 
14,2000) 

C TRANSFERS  INFO  OF  A GIVEN  PRISM  INTO  WORKING  ARRAYS  AND  FORMS  INITIAL  GUESS 
C TO  AID. 

C IMF  NS  I ON  X ( NUMONC  ) , Y ( NUMOND ) , B ( MUA ) , VALS t NOEQ , NLMOND ) 

C DETERMINE  NUMBER  OF  NODES  OF  PRISM — NDPPRM 
NCPPRM=1 
MK  = NUMONC— I 
CO  20  J=i,MK 

IFINRIJ, NONUM). EQ.0)G0T0  25 
20  NCPPRM  = NCPPRM-»1 

C TRANSFER  CATA  INTO  X.Y.VALS  FROM  XY.VAR  AND  SET  TO  ZERO  UNUSED  PORTION  OF 
C XfYtVALS 

C CENTER  NODE  INFO  IN  COORDINATE  SYSTEM  WITH  ORIGIN  AT  CENTER  NODE 
25  X(1I=0.0 

vm-o.o 

CO  35  1=1, NOEQ 

C KOLC  REFERS  TO  INFO  AT  KNOWN  INITIAL  LEVEL 
35  VALSI I»1)*VARIK OLD, I, NONUM) 

C NEIGHBORING  NODE  INFO  IN  COORDINATE  SYSTEM  WITH  ORIGIN  AT  CENTER  NODE 
CO  30  K = 2 ,NCPPRM 
JSUB=NRIK-1, NONUM  ) 

X(K)=XY! 1,JSUB)-XY( 1, NONUM) 

Y(K)=XYI2,JSUB)-XY(2, NONUM) 

CO  30  1=1, NOEQ 

30  VALSI  I , K l=VAR(KOLC,I,JSUB) 

C IF  NUMBER  OF  TRIANGLES  IS  LESS  THAN  MAX  NUMBER  ALLOWED  IN  A PRISM,  SET 
C RFMAINCER  OF  X.Y.VALS  TO  ZERO 

IFINCPPRM.EQ.NUMONOIGOTO  40 
JIM=NCPPRM*1 
CO  45  K=JIM, NUMONC 
X I K ) =0 . 0 
Y ( K ) =0.0 
CO  45  1=1, NOEQ 
45  VALSI  I,K  )=C.C 

C ASSUMPTION  MUA/NOEQ  IS  AN  INTEGER  I PARAMETERS  A ARE  EQUALLY  SPACED  AMONG 
C UNKOmNS) 

40  MCIV=MUA/NOEC 

C FORM  INITIAL  GUESS  OF  PARAMETERS  All) 

C IN  THE  TRANSLATEC  COORDINATE  SYSTEM.  WE  TAKE  A 1 1 , 2 , 4 ,5 , 7 , 8, 10  AND  11)=0 
C AND  Al I*3)=VARTI I, NONUM) 

MCIVM1  = MC  IV- 1 


S7 


oooo  ooooo  n n n r>  n 


CO  50  I=1»N0EQ 
P < HO  I V* I >«VARTI l.NONUM) 

CO  50  J»1,MCIVMI 
50  BINCIV*I-J  1=0.0 
RETURN 
ENC 

SUBROUTINE  SOL VER ( NUMOND , NOEQ ,MU A. T I ME ,DT , NDPPRM ,NONUM, X , V , VALS , B , 

IFU, FUSPEC) 

THIS  SOLVER  IS  THE  BASIC  NEWTON-RAPHSON  SOLVER 
COMPUTES  THE  FINAL  VALUE  OF  THE  PARAMETERS  All)  BY  THE  NEWTON-RAPHSON  METHOD 
WITH  THE  TERMINATION  CRITERIA  GIVEN  BY  THE  SUBROUTINE  TERMIN 
FU  ANC  FUSPEC  ARE  SUBROUTINES  TO  COMPUTE  NECESSARY  QUANTITIES  FOR  NEWTON- 
RAPHSON  ITERATION  PROCEDURE 
EXTERNAL  FU, FUSPEC 
LOGICAL  INDICT, PROCED 

DIMENSION  B(MUA) , X ( NUMOND  I , Y ( NUMOND ), VALS ( NOEQ , NUMOND ) 

Cl  MENS  ION  CAJVNII 12,12),F( 12), SAVE  I 12), A (12) ,WTS(12) ,DMCRFA(12) 

HXITPJ  IS  THE  MAX  NUMBER  OF  ITERATIONS  ALLOWED  PER  JACOBIAN  EVALUATION 
MXNOJE  IS  THE  MAX  NUMBER  OF  JACOBIAN  EVALUATIONS  ALLOWED. 

IF  PR0C6C=T , PROGRAM  WILL  CONTINUE  EVEN  THOUGH  CONVERGENCE  CRITERIA  IS  NOT 
MET  AFTER  ALL  ALLOWABLE  ITERATIONS. 

IF  PROCE  C = F , PROGRAM  WILL  STOP. 

CATA  MXITPJ.MXNOJE/l, 10/,  PROCED/.F AL SE . / 

CO  20  1=1, MUA 
SAVE!  I ) = P(I) 

20  A<  I»  = B( I I 

THE  NEWTON-RAPHSON  METHOD 

***THE  MAX  NORM  IS  USEC  FOR  THE  CONVERGENCE  CRITERIA  SINCE  AID  ARE 

C I MENS  IONAL , CMCRFA(I)  IS  USEC  TO  MAKE  THEM  NON-DIMENSIONAL  FOR  TESTING 
CONVERGENCE 

CALL  SC ALFAINUMONC, NOEQ, MUA, DT, VALS, DMCRFA) 

CO  30  L=l, MXNOJE 

C COMPUTE  THE  VECTOR  F AND  JACOBIAN 

CALL  FU( NUMOND, NOEQ, MUA, T IME , CT , NDPPRM , X, Y , V AL S , A , F ,C A J VN I ,RES) 

C COMPUTE  THE  INVERSE  OF  JAC-CAJVNI 
C MAT INV  IS  SUBPROGRAM  AVAILABLE  ON  BRLESC  SYSTEM 
CALL  MATINV(CAJVNI,MUA,F,MUA,C,DET) 

IF  (CET  .EQ.  0.0)  GO  TO  35 

C ITERATE  MXITPJ  TIMES  BEFORE  EVALUATING  JACOBIAN  AGAIN 
CO  30  K=i, MXITPJ 

C CALCULATE  NEW  VALUES  OF  PARAMETERS 
CO  45  1=1, MUA 
CO  45  J= 1 , MU A 

45  At  I l>A(1  l-CAJVNII  I,J  )*F( J ) 

C CHECK  FOR  CONVERGENCE 

C STPCRT ,RELERR,INCICT,SAVE,WTS  ARE  OEFINED  IN  TERMIN 

CALL  TERMIN  I NUMONC, NOEQ , MU A , CT .NDPPRM , INDIC T , STPCRT .RELERR , VALS , 

1 A, SAVE.WTS, DMCRFA, NONUM) 

C IF  AID'S  ARE  WITHIN  TOLERANCES,  INDICT*TRUF 
IFIINCICT)  GO  TO  50 
C ITERATE  AGAIN 

CO  55  1=1, MUA 
55  SAVE!  I ) = A I I ) 

C ♦♦FUSPEC  COMPUTES  ONLY  F OF  NEWTON-RAPHSON,  NOT  THE  JACOBIAN 

IF  IK. LT. MXITPJ)  CALL  FUSPEC I NUMOND, NOEQ, MUA , T I ME  ,0T .NDPPRM , X , Y , VAL 
IS, A, F, RES) 

30  CONTINUE 
C NO  CONVERGENCE 

NOITER“MXlTPJ*MXNOJE 
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PRINT  HO 

PRINT  81.N0NUM, TIME,DT,NOITER,STPCRT,RELERR,PROCED 
PRINT  82 

PRINT  83,(B( I),A(  I I,  UTS  (I).  I >1,  HUAI 
IFI .NOT.PROCECISTOP 
C CONVERGENCE  IS  OETAINEO 
50  CO  HO  1=1, MUA 
60  ei n*Ai 1 1 

PRINT  92.N0NUM.L 

92  FORMAT! 1H , 15X , 'NOCE • . I 5, 3X . 1 4 , • I TERAT I ONS • I 
RETURN 

C TROUBLE  WITH  FINCING  THE  INVERSE  JACOBIAN 
35  PRINT  91 

PRINT  85,  NON UH, TIME, DT.L 
PRINT  86 

PRINT  87,  (XI I)  ,Y( I 1 , I=l,NUMONDI 
PRINT  90 
PRINT  88 

CO  47  1=1 , NUMOND 

47  PRINT  89,I,(VALS(J,I  ),J=1,N0EQ) 

STOP 

80  FORMAT! 1H,//,  • NO  CONVERGENCE  • I 

81  FORMAT ( 1H, • NUM8ER  OF  NOOE  IS',15,'  T I ME= • ,F 12 . 7 , • DELTA  T=',F8.5 

l,//,'  NUMeER  OF  ITERATIONS** , 14, • CRITERION  FOR  CONVERGENCE  IS'. El 
24.5,//  • LAST  COMPUTED  RELATIVE  ERROR  IS*,E14.5,  / ,*  PROGRAM 

3WILL  NOT  STOP— *,L3) 

82  FORMAT  (IF,*  INITIAL  A(II  VALUES  * , 5X , *L A ST  All!  VALUES', 7X, 

1 • HE l GHTS ' I 

83  FORMAT  (1H,  3X , E 14 . 5 , 8X , E 14 . 5 , 4X,E 14 . 5 I 

85  F0RMATI1H,*  NUMBER  OF  NODE  IS',15,*  T IME* • ,E 15 . 8 , • DEL  T=',E15.8 
l,/'  JACOBIAN  WAS  eEING  EVALUATED  THE',I3,'TH  TIME'I 

86  FORMAT  (IF,  3X ,' X- VALUES ', 1CX ,' V- VALUES • I 

87  FORMAT  (IF,  E 14 . 5 , 4X , E 14 . 5 ) 

88  FORMAT  (IF,*  NODE  NUMBER • , 3X , • U- VEL OC I T V • , 3X , • V- VE LOC I T Y • ,4X , • PRES 
1SURE' ,4X, 'INT. -ENERGY • ,3X, 'SWIRL  VELOCITY* I 

89  FORMAT  (IF,  5 X , 1 3 , 6X , E 12 . 5 , IX , E 12 . 5 , IX ,E12 . 5 , 2 X , E 12 . 5 ,2 X ,E 12 . 5 1 

90  FORMAT  (IF,'  VALUES  AT  INITIAL  SURFACE  (0T*0)  ARE  'I 

91  FORMAT  (IF,  //,'  STOP  IS  CUE  TO  DET*0 .0  IN  M AT  I N V-SOL VER • I 
ENC 

SUBROUTINE  T E RM IN ( NUMOND, NCEQ, MU A, DT.NDPPRM, INDICT, ST  PORT .RELERR, 

VALS, A,SAVE,WTS,DMCRFA,NCNLM) 

C USES  MAX  NORM  OVER  I=1,...,MUA  TO  DETERMINE  CONVERGENCE 
LOGICAL  INC1CT 

C IMF  NS  ION  VALS(NOEQ,NUMOND),A(MUA),SAVE(MUAI ,WTS(MUA) 

C I MENS  ION  CMCRFA(  1 2 ) , EX  T ( 12) 

C V>  T S I 1 I IS  THE  WEIGHT  ASSOCIATED  WITH  AID 

C WTS  ALLOWS  MODIFICATION  OF  MAX  NORM  TO  A WEIGHTED  MAX  NORM 
C STPCRT  IS  UPPER  BOUND  OF  ALLOWABLE  RELATIVE  ERROR 

C THIS  PROGRAM  FAS  A DATA  STATEMENT.  CHANGE  DATA  TO  AGREE  WITH  DIMENSION 
C OF  ARRAY  WTS(MUA). 

DATA  EXTSTC/l.E-5/,(EXT( I), 1=1, 121 /0.0,0.0,l.0,C. 0,0. 0,1. 0,0.0, 
10.0, 1.0, 0.0, 0.0, 1.0/ 

STPCRT=EXTSTC 
CO  20  1=1, MUA 
20  WT S ( I l = EXT ( I I 

C COMPUTE  THF  MAX  RELATIVE  ERROR 

RELcRR=(CMCRFA(  II*  WTS( 1 > * I A ( 1 1-SAVE ( 1 1 ) 1**2 

CO  40  1=2, MUA 

S T= ( CMCRF A ( I ) *WTS ( I ) * ( A I I )-SA VE ( I I ) )**2 
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C KFLERR’RELAT  IVE  ERROR 

40  RELFRR*MAX1F(ST,RELERR) 

PELERR*SQRTF(RELERR) 

C COMPARISON 

C INDICT  HAS  VALUE  TRUE  IF  CONVERGENCE,  ELSE  FALSE 
INC!CT=RELERR.LT.STPCRT 
RETURN 
ENC 

SUBROUTINE  SC ALFA  I NU MONO , NOE 0 ,MUA,OT , VALS, SF ) 

COMPUTES  DIMENSIONAL  CORRECTION  FACTORS  FOR  EACH  AID  SQUARED 

FOR  EXAMPLE  U-U0+ ( A ( 1 ) *X* A ( 2 > *Y+A ( 3 > ) *T  UNITS  FOR  A(1),A(2>  AND 
A ( 3 ) ARE  1/ ( SEC*SEC  ) , l/(SEC*SECI,  ME TRE 5/ 1 SEC* SEC ) SO  THE 
EVERY  TERM  HAS  UNIT  (METRES  PER  SECOND) 

Cl  MENS  ION  V AL  S ( NOEQ, NUMONC ) , SF ( MU A ) 

CHARACTERISTIC  VELOCITY  IS  SOUNC  SPEED  AT  DT*0  ANO  CENTER  NODE 
CALL  SNCSPDIVALSI 3,1  ), VALS (4,1) ,REFVEL> 

CHARACTERISTIC  LENGTH  IS  0T*SNDSPC 
REFLNG=CT*REFVEL 

CHARACTERISTIC  PRESSURE  AND  INTERNAL  ENERGY  ARE  THOSE  VALUES  AT  CENTER 
NOCE  ANC  CT*0 

REFPR  = VALS( 3,  l ) 

REFIE*"ALS(4,1) 

COMPUTE  THl  CIMENSIONAL  CORRECTION  FACTORS  FOR  THE  AID'S 
SF ( 1 1 *DT  *CT 
S F ( 2 ) = SF ( 1 ) 

SF ( 3 ) * CT/REFVEL 
SF(4)*SF(1) 

SF(5)*SF(2) 

SF(6)=SF( 3) 

F=  CT  *REFLNG 
SF ( 7 ) *F/REFPR 
SF ( fl ) =SF ( 7 I 
SF ( 9 ) *DT  / PE F PR 
SF ( 10  )=F/R£F I E 
SF ( 1 1 ) = SF ( 10  ) 

SF ( 1 2 ) =CT  /REF  IE 
IF(N0EQ.NE.5 (GOTO  30 
SF ( 13  ) = SF ( 1 ) 

SF ( 14  ) *SF ( 2 ) 

SF ( 15  l = SF ( 3 ) 

30  RETURN 
ENC 

SUBROUTINE  FU ( NUMONO ,N0EQ  , MUA , T IME ,0T ,NOPPRM , X , Y , V ALS , B ,F , F A ,RES ) 

COMPUTES  TRIPLE  INTEGRAL ( F ( I I , I *1-MUA ) OVER  A PRISM  CF  THE  SUM  OF  TERMS  0 ( K ) * 
D ( K , I ) ,K* l-NO EQ  ANC  ITS  FIRST  PARTIAL  DERIVATIVES  (FA(I,J)),IF  DESIREO 
D ( K I IS  THE  RESICUAL  CORRESPONDING  TO  THE  KTH  EQUATION 

C IMF  NS  I ON  X ( NUMONC ) ,Y(NUMOND) ,B( MUA ),VALS( NOEQ.NUMOND) ,F( MUA), 
IFA(MUA.HUA) 

Cl  MENS  ION  RHC(9, 10),PF( 12),PFA( 12,121 

IN  SOLVING  A ( NEW  ITERATE  ) *A < OLD  I TER ATE ) - 1 NVERSE ( J ACOBI AN ) *F , UPDATE 
JACOBIAN  ANC  F IF  INDI*0,  ONLY  F IF  INDD1 
CO  20  1*1, MUA 
CO  20  J* 1 »MU A 
20  FA(I,J)*0.0 
INC  1*0 

C***«  RES  IS  THE  LEAST  SQUARE  RESIDUAL  THAT  IS  TO  BE  M IN  AMI  ZED 
RE  S = 0 .0 
GOTO  25 
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ENTRY  FUSPEC ( NUMOND, NOEQ, MUA, TIME.DT.NOPPRM, X , Y , VAL S ,B ,F ,RE  S I 

1 NC  I - 1 

RES=0.0 

25  CO  30  I * 1 , MU A 
30  F(  11*0.0 

C APPROXIMATE  TIME  INTEGRAL  BY  A TWO  PT.  GAUSSIAN  QUADRATURE 
JUDGE  *0 

C FIRST  GAUSSIAN  TIME  POINT 
GAUST*0.21132A865A*DT 

C CALCULATE  VARIABLE  FROM  EQUATION  OF  STATE  AND  ITS  DERIVATIVES  AT  PRISM  NODES 
35  CALL  ARECST  ( NUMONO ,NOEQ , MUA , NDPPRM ,GAUST , X , Y , VAL S , B ,RHO I 
C CALCULATE  SPATIAL  INTEGRALS  OVER  PRISM  FOR  A GIVEN  GAUSSIAN  TIME 
CO  AO  K*  I iNOEQ 

C * * * * PF  ANC  PFA  CORRESPOND  TO  THE  SPATIAL  INTEGRALS  ASSOCIATED  WITH  F ANC  FA  AT 
C A GIVEN  GAUSSIAN  TIME 

CALL  EVALPR( NUMOND, NOEQ, MUA .NDPPRM, IND I , GAUST , X , X , Y , VAL S , B , RHO.PF , 
1PFA.PRES I 
RES*RES*PRES 
CO  AO  L * 1 » MU  A 
F(L)*F(LMPF(L> 

C IF  I N C I = I DO  NOT  CALCULATE  FA 
IF  ( INCI  .EC.  1 ) GOTO  AO 
CO  AO  j*;,mua 
FAIL#  Jl*FA(L,  J I»PFA(L,  J» 

AO  CONTINUE 

C SFCONC  ANC  LAST  GAUSSIAN  POINT 
IF  (JUDGE  .NE.  01  GOTO  A5 
JUCGE*  l 

GALST«0.78e67513A6*CT 
GOTO  35 

C MULTIPLY  SUMS  BY  CONSTANT  HEIGHT  FACTOR  0.5*DT 
RES=0.5*CT*RES 
A5  CO  50  L * 1 » MU  A 
F(L)=0.5*CT*F(L  ) 

C IF  INCI  = I L'O  NOT  CALCULATE  FA 
IF  (INCI  .EC.  11  GOTO  50 
CO  55  J*  1 , MU  A 
55  FA(L,JI*0.5*CT*FA(L,J) 

50  CONTINUE 
RETURN 
ENC 

SUBROUTINE  E VAL PR ( NUMOND , NOE Q ,MUA , NDPPRM , I ND I ,GALST,K,X,Y,VALS,B, 

1 RHO.PF, PFB.PRESI 

C COMPUTES  SPATIAL  INTEGRALS  OVER  A PRISM  FOR  A GIVEN  K RESIDUAL  TERM  AND 
C GAUSSIAN  TIME 

COMMON/ INARRS/XTRI ( 3I.YTRI ( 31 ,U(3) ,V( 3) ,P( 31 ,E( 3) ,SV(3) .RHOTRI (3,1 
1 0 ) /OUT ARKS /D( 3),0A(  12,3)  ,DAA ( 12 , 12 , 3 I 

C IN  COMMON/OUTARRS/.WE  HAVE  MUA*12.  WHEN  MUA  IS  NCT  12,  MANUALLY  CHANGE 
C DIMENSIONS  IN  OUTARRS 

Cl  MENS  ION  X ( NUMONC  I, V( NUMOND  ) , VAL S ( NOEQ, NUMOND ) ,RHO( NUMOND, 101 ,B(M 
1 UAI,PF(MUAI,PFB(MUA,MUA),RR(3I,RR2(3)  , RR3 ( 3 I , RRA ( 3 ) 

C SV  DENOTES  SWIRL  VELOCITY  (N0EQ=5> 

PRES*0.0 
CO  10  I * l ,MU A 
10  PF (11=0.0 

C IF  INCI=1,THEN  JACOBIAN  IS  NOT  TO  BE  CALCULATED 
I F ( I ND I .EC.  1)  GO  TO  15 
CO  16  1*1, MUA 
CO  16  J*  1 ,MUA 
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16  p f e 1 1, j)*o.o 

C COLLECT  INFO  FOR  * GIVEN  TRIANGLE  INTO  XTRI , YTRl ,U , V ,P ,E , SV .RHOTRI 
15  XTRim=>X<ll 
YTRim*Y(il 
UIlt«VALS(l,n 
vm*vALS(2,n 
P ( 1 ) 1 V AL  S ( 3 » 1 ) 

E( 1)»VALS(A,1) 

SVI1I>0.0 

IF (NOEO.EC.5 ) SV( 1 )*VALS( 5 t 1 ) 

NUMTRl=NCPPRM-l 

C LOOP  OVER  TRIANGLES  COMPOSING  THE  PRISM 
CO  20  ITH  I* l.NUMTR I 

C COLLECT  REMAINING  CATA  FOR  A GIVEN  TRIANGLE 
l V2= I TR I ♦ 1 
I V 3=  ITRl*2 

IFIITRI.EC.NUMTRI  )IV3*2 
XTR  I ( 2 I *X ( IV2  ) 

XTRII3I-XI IV3I 
YTR  l ( 2 ) = Y ( IV2  ) 

YTR I I 3 I = Y I I V 3 I 
UI2I-VALSU.  IV2) 

U<3)=*VALS<1,IV3) 

V(2)*VALS(2.IV2) 

V(3)=*VALS(2,IV3I 
P(21*VALS(3f IV2) 

PI 3)=VALS( 3, IV3  > 

E I 2 I *VAL  S I A. IV2I 
E(3)=VALS(A, IV3 I 
IFIN0EQ.EC.5 IGOTO  25 
S V 1 2 ) *0 . 0 
S V I 3 I *0 .0 
GOTO  26 

25  SVI21*VALS(5,  IV2) 

SV(3)=VACSI5,IV3) 

26  CO  27  J=l,10 
RHOTRI  C 1, J)*RHO< 1, J I 
RH0TRII2, J)=RHO< IV2,J) 

27  RHOTRII 3, J l*RHO( IV3, J) 

C COMPUTE  QUANTITY  PROPORTIONAL  TO  TRIANGLES  AREA 
CET  = CETERM(XTRI,YTRI  ) 

C CSUB  COMPUTES  KTF  RESIDUAL  ID)  AND  ITS  DERIVATIVES  AT  VERTICES  OF  GIVEN 

C TRIANGLE  (CA,CAA) 

CALL  CSUB(MUA,GAUST, INOI ,K,OET,B I 

C PROCEEC  TO  CALCULATE  THE  SPATIAL  INTEGRALS 
CALL  AREAINIC.CtOET, SOROS) 

PRES*PRES«SOROS 
CO  30  I * 1 1 MU A 
CO  35  J=l,3 
35  RR(J)*DA<  l,J) 

CALL  AREAINIC.RRfDET.SXY) 

30  PFI I l«PF( II+SXY 

C IF  IN0I=1  THEN  JACOBIAN  IS  NOT  TO  BE  CALCULATED  INCREMENT  ITRI. 

IF!  INDI.EQ.l IGOTO  20 
CO  AO  1 = 1, MUA 
CO  A2  L*i,3 
A 2 RRA | L 1 = 0 A ( 1 1 L I 
CO  AO  J=I,MUA 
CO  A 7 L*  1 , 3 
RR2IL  l=OA(J,L  I 
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47  RR3(LI=CAA(  I , J , L I 

CALL  AREA  I N ( C,RR3, GET, SXY1 ) 

CALL  AREAIN(RR4,RR2,0ET,SXY2I 
40  PFe( I,J1=PFB(  I,J»«SXY1>SXY2 
20  CONT INUE 

C REFLECT  CALCULATED  UPPER  TRIANGULAR  PORTION  OF  PFB  TC  OBTAIN  PFB.IF  NECESSARY 
IF  ( INCI  .EG.  1 1 GO  TO  101 
PMUA*MUA-1 
CO  SO  l-l.PMUA 
I I 1 = 1*1 

CO  SO  J* I 1 1 » HU A 
SO  P F E ( J * I I =PFB ( I, Jl 
101  RETURN 
ENC 

SUBROUTINE  CSUBIMUA.T, IND I , K SPEC ,DET , A J 
C I PENS  ION  A ( 1 2 ) 

C IPENS ION  WTSI4) 

COPMON/INARRS/AX<  3 I , AY ( 3 I , AU ( 3 I , AV ( 3 I , AP ( 3 I , AE ( 3 I , ASV ( 3 I , ARHOI 3, 10 
1 l/0UTARRS/0( 3>,CA( 12 , 3 I . DA A ( 1 2 , 12 , 3 I 
C WTS  ALLOW  FOR  A WEIGHTED  LEAST  SQUARES  RESIDUAL  FORPLLATI ON  IF  DESIREO 
CATA  (WTS<  I )• 1 = 1*4)  /l. 0,1. 0,1. 0,1.0/ 

C nsue  COMPUTES  THE  VALUE  OF  THE  RESIDUALS  O AND  ITS  DERIVATIVES  (DA.DAAI  WRT 
C a ( j ) AT  THE  VERTICES  OF  THE  TRIANGLE  BEING  CONSIDERED. 

C TC  MAKE  THE  INDIVIDUAL  TERMS  CF  THE  LEAST  SQUARES  RESIOUALS  NON-DIMENSIONAL 
C CIVICE  FACH  TERM  BY  AN  APPROPRIATE  FAC TOR-DMCF  SEE  THE  END  OF  EACH  SUBSECTION 
C KSPEC=1,2,3,4 

C rsue  IS  CODEC  PRESENTLY  FOR  MUA*12  AND  K*l,2,3,  OR  4.  DSUB  CHECKS  FOR  THESE 
C VALUES 

IFIPUA.NE.12IG0T0  101 

C ON  EACH  TRIANGLE  USE  A LINEAR  APPROXIMATION  TO  U.V.P.I  AT  TIME  LEVEL  CORRES- 
C PONCING  TO  DT*0.  BELOW  ARE  FACTORS  NEEOED  FOR  THIS  APPROX. 

FX 1 = AX ( 3 )-AX( 2 IS  FX2=AX(3I-AX(IIS  F X 3 = AX  I 2 ) - AX C l I 

FYl=AY(3l-AY(21S  F Y2* A Y ( 3 1- A Y ( 1 I $ FY3=AY( 2 1-AY (11 

IF  IKSPEC  .EC.  II  GO  TO  15 

IF  ( KSPEC  .EC.  2)  GO  TO  30 

IF  (KSPEC  .EC.  31  GO  TO  45 

IF  (KSPEC  .EC.  41  GO  TO  60 

PRINT  99, KSPEC 

STOP 

C PEGIN  CALCULATIONS  FOR  K=1  (CONTINUITY  EQN. I 

15  UOX=(  ( AU ( 2 l-AU( 3)  I • AY ( 1 l-(AU( 1 I-AU 1 3 1 I *A Y ( 2 ) ♦ ( AL ( 1 1- AU (2  I I *AY (3)1 
1 /CET 

EOX= ( ( AE(2)-AE( 3) )*AY( 1 I- ( AE ( 1 1-AE ( 3 1 ) *AY < 2 1 ♦ ( AE ( 1 1- AE < 2 I I • AY ( 3 II 

1 /CET 

POX=(  (AP(2  I — AP ( 31  |*AY( 1 l-(AP( 1 l-AP ( 3 1 I * A Y ( 2 I ♦ < A P ( 1 1 - AP ( 2 I I *AY ( 3 I I 
1 /CET 

VOY=(-( AV(2I-AV(3»  I * AX ( 1 1 ♦ ( A V ( 1 I- A V ( 3 1 1 * AX ( 2 I - ( A V ( 1 1 - AV ( 2 I I * AX ( 3 1 I 
1 /CET 

E0Y=(-(AE(2I-AE(3I l*AX( 1 1 A ( AE ( 1 l-AE I 3 I I MAX  1 2 I- ( AE  < 1 l-AE (2  I 1* AX (31  I 
1 /CET 

P0Y=(-(AP(2»-AP(3l  >*AX(1 !♦( API  1 l-AP ( 3 II *AX I 2 I- ( AP ( 1 1 -AP ( 2 I »*AX(3> I 
l /CET 

PX»P0X*T*A(7I*  PY  *P0 Y *T*A ( 8 I 
EX»E0X*T*A(10I*  EY»E0Y4T*A( 111 
CR«U0X4V0Y4T*I  A(  1HA(  51  I 
TTT=T*T 
CO  20  1*1,3 
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X- AX  I I ) $ Y-AYI  1 I 

R-ARH0II.1)*  RP-ARHOI 1,2)*  R I ■ ARHO ( I « 3 I 
RFP-ARHOI I, A It  RII-ARHOI 1,5)*  RPI -ARHOI I ,6 1 

RPPP= ARHO ( I » 7 ) $ RlII-ARHOI 1,0)*  RPPI -ARHOI I , 9 ) t RPI I -ARHO 1 1 , 10 > 

UT*»(ll»X*AI2)mA(3l 

VT*A(AI*X*A(5I*Y+A(6) 

U-AUI I I ♦ T *UT 
V-AVI  1 1 ♦T*VT 
PT-A(7)*X*A|8)*Y*A|9) 

ET*  A ( 10 ) *Xt A I l 1 )*Y*A I 12) 

CRP=PT*l)*PX*V*PY 

CRI=ET*U*EX>V*EY 

C TFF  FIRST  RESIOUAL  D AT  THE  VERTICES  OF  THE  TRIANGLE  <I»1,2,3IIS 
Cl I)-RP*CRP*Rl*CRI*R*CR 

C IFE  DERIVATIVE  OF  THE  1ST  RESIOUAL  WRT  AIK)  IS  OAIK.I)  AT  THE  VERTICES 
C (1=1,2,31 

RTT  = R*T 
XPTTU-X*  T *U 
VPTTV-Y*T*V 
STOR=IRP*PX»Rl*EX)*T 

CA(l,I)=STOR*X«RTTt  OA I 2 , I ) -STOR*Yt  OA(3,II«STOR 
STOR«(RP*PY  + R I*EY  l*T 

CAIA,  I l»STOR*Xt  DAIS, I ) * STOR * Y + k TT t DA|6,II-ST0R 
STOR=(RPP*CRP«RPI*CRI*RP*CR l*T 

CA( 7, I )-RP*XPTTU+STOR*Xt  0A( 8,1 ) «RP* YP TTV*S TOR* Y * OA (9, I I »RP*STOR 
STCR=(RPI*CRP*RII*CRI*RI*CR)*T 

CA(10,I I=RI*XPTTLMSTOR*XtOA(ll,  I )-RI*YPTTV*ST0R*YtDA(12»I )*Rl«-STOR 

C 2NC  CERIVATIVES  OF  0 (OAAIL.K.I))  ARE  NOT  TO  BE  CONPLTEO  IF  INDI-l 
IF! INCI.EC.l  I GO  TO  20 

OAAIL.K.I)  IS  THE  DERIVATIVE  OF  THE  1ST  RESIOUAL  0 WRT  AIL)  AND  AIK)  AT  THE 
VERTICES  (1=1.2.31 
XTT=X*Tt  XTT 2-XTT *T 
YT  T = Y*T  t YTT2-YTT*T 
RPTT2-RP*TTT*  R ITT2-R  I*TTT 
XPTUTT=XPTTU*Tt  YPTVTT  -YPTTV*T 
SPITXU=RPI*XPTUTTI  SP  I TY V-RP  I * YPTVTT 
GPXT-(RPP*PX«RPI*EX)*TTTt  G I XT  * I RP I *PX*R 1 1 *E X I *T TT 
GPYT-(RPP*PY*RPI*EY)*TTTt  GIYT-(RPI*PY*RII*EY)*TTT 
CO  2 A L-1,6 
CO  2 A H-L  * 6 
2A  CAAIL.M, I » - 0 . 0 

ST0RA»GPXT*X«RPTT2 
ST0RB-GIXT*XfRITT2 

CAAI 1,7, I )»( ST0RA«RPTT2l*Xt  DA A 1 1 , 8 , I > -STOR A*Y * CAAI 1 ,9, 1 l-STORA 
CAAI 1,10, 1 I- 1 ST0Re*RITT2 »*X*  CAAI 1,11,1 ) = STORB* Y *OAA ( l , 12 , I ) -STORB 
STOR  = GPXT  *Y 

CAAI2.7, I )«STORA*Yt  DA A I 2 , 8 , I I * STOR • Y$  OAA 1 2 ,9 , 1 »« STOR 
ST0R-GIXT*Y 

CAAI 2, 10, I)-STORB*Y$  OAA I 2 , 1 1 , I ) -STOR* Yt  OAA (2 , 12 . 1 ) -STOR 
CAAI 3,7, I )«STORAt  DA Al 3, 8 , I ) -GPXT*Yt  OAA I 3, 9 , I) -GPXT 
CAAI 3,10, I l-STORBI  0 A A I 3 • 1 1 , I l-STOR I OAA I 3, 1 2 , I ) «G I XT 
ST0RA-GPYT*Y>RPTT2 
ST0RB-GIYT*Y*RITT2 
S TOR-GPYT  *X 

C A A I A , 7, I )-STOR*XS  DAAI A, 8. I )-STORA*Xt  OAA I A ,9 , I I -STOR 
S TOR-G  I YT *X 

CAA(A,10, I l-STOR*Xt  CAAI A, 11, I )* STORB* Xt  OAA I A , 12 , 1 1 -STOR 
CAAI 5,7, I ) - STOR A*X  $ 0 A A I 5 , 8 , I ) * I STOR A*RPTT2 ) *Y t DA A 1 5 ,9 , 1 ) -STOR A 
CAAI 5,10, I)-STORB*Xf  OAA I 5, 11, I )>( STORB*R I TT2 ) * YIOA A 1 5, 12 , 1 ) -STORB 
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CAA16.7,  I >-GPYT*X»  OAA ( 6 , 8 , l ) -STORAA  DAA (6 ,9 , 1 ) -GP VT 

C A A ( 6 , 1 0 , I ) » S TOR  t DAA(6, 11, II-STORBS  DAA < 6 , 1 2 , II «G I YT 

STORC*RPPP*CRP«RPPl*CRlARPP*CR 

STORO«RPPI*CRP*RPI l*CRI*RPI*CR 

STOR A»RPP*XPTUTT 

STORB-STORCAXTT2 

CAA( 7,7, I I«(2.*STCRAaST0RBI*X 

CAA<  7,8, 1 l»(STORA«STORB)*Y'»RPP*YPTVTT*X 

C A A ( 7 ,9, I |*STORA^STORB*RPPAXTT 

STORA-SPITXU 

STORB«STORC*XTT2 

CAA( 7,10, n»(2.*ST0RA*ST0RB  »*X 

CAA( 7,11, H»( STORA  + STORB  >*VASPITYV*X 

CAA( 7, 12, I)*STORA«STORBARPIAXTT 

STORA«RPP*YPTVTT 

ST0RB»ST0RC*YTT2 

CAAC8.8, I )*(2.*ST0RA*ST0RB)*Y 

CAA( 8,9,  I >«STORA*STORBARPP*YTT 

STORA*RP!*YPTVTT 

ST0RB*ST0RC*YTT2 

CAA( 8,10, 1 )«(STORA*STORBI*X*SPITXU*Y 

CAA(8,ll,  l»»( 2.AST0RA*ST0RB1*Y 

C A A ( 8 , 1 2 , I )*STORA«STOR0*RPI*YTT 

CAA(9,9,  l »*(2.0ARPP*ST0RC.AT  l*T 

STORA*(STCRC*T  + RPI  )*T 

C A A ( 9 , 10 , I )=STORA*X*SPITXU 

CAA(9,11, I )=STORAAY*SPITYV 

CAA ( 9, 12,  I (*STORAARPIAT 

STORE-RPI  I*CRP*RI  I I *CR  I *R  1 1 *CR 

S TOR  A*R I MXPTUTT 

ST0RB*ST0REAXTT2 

CAA( 10, 10,I)=(2.OASTORA*ST0RB>*X 

CA Al 10,11 , I )‘l STORAaSTORB )*V  + X*R1 I*YPTVTT 

C A A ( 10,12,1 )*STORA*STORB*RII*XTT 

STOKA=YPTVTT*RI  I 

STCRB»STORE*TTT 

CAA (11, l 1,1 >=(2.*ST0RA+STCRBAY)*Y 
CAA I 11, 12, I > = ST0RA»ST0RB*YAR1  I*YTT 
CAA(12,12,n  = 2.0*RlI*TAST0RB 
CO  28  M= i , 1 1 
K P l = K ♦ 1 

CO  28  L=KP1,12 
28  CAAIL.H, I I»CAA(P,L,  I ) 

20  CONTINUE 

C FOR  CONTINUITY  ECUATICN  THE  CHARACTERISTIC  QUANT  ITY (CPCF I IS  OENSITY  AT 
C T*C  ANO  CENTER  NODE 

CHCF*ARHO<  1,  1 I 

co  200  ic-i,a 

C(  IC l*WTS< 1 >*C<  IC  I7CMCF 
CO  200  JC>1,RUA 

CA  CJC,  IC )«WTS< 1 > *0A( JC,  ICI/OPCF 
CO  200  KC*1,HUA 

200  CAAIKC,  JC,  IC  )«HTSm*DAA(KC,JC,  IC  I/OMCF 
RETURN 

C CECIN  CALCULATIONS  FOR  K-2  IX-NOHENTUH  EQN.) 

TO  UOX=(  ( AU(2)-AU(3I )» AY  III- 1 AO « 1 l-AU <3 M *AV 1 2 I ♦ I AM  1 1 > AU (2 > I » AY (31  I 
1 /CET 

UOV>  l-(AU(2l-AU(3M*AX(  1>«(  AU(  1 I-AU  (3  M *AX  I 2 I- I AL  ( 1 1- AU  C2  I I *AX  1 3 II 
l 7CET 
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P0X=1  (AP<2)-AP(3))*AY<l)-(AP<l)-AP(3))*AV(2IMAP(l)-API2>)*AY(3>> 

1 /CET 

UY *UOY*A  I 2 > *T  $ UX»UOX*AI l)*T 

UYTTT*UY*T*T 

CR 1*  1 »0«T  *UX 

CR  1TT*CR  i*T 

CO  35  1*1,3 

X * AX  I It  t V*  AY  ( I I 

R* ARHOI I , 1 ) S PP * ARHO ( 1 , 2 ) t R I * ARHOI 1,3) 

PFP=ARHO( I, Alt  RII*ARHOI I,5)$  RPI*ARHO( 1,6) 

UT*A(I)*X*AI2)*Y4AI3> 

VT*A(4)*X«A(5)*Y«AI6) 

U*  AIM  1 I ♦TAUT 
V-AVI  I)*r*VT 
CR*UT»U*UX*V*UY 

C THE  SECOND  RESIOUAL  0 AT  THE  VERTICES  OF  THE  TRIANGLE  11*1,2,3)  IS 
Cl  I)*R*CR«POX4AC  7 )*T 

C THE  CERIVATIVE  OF  THE  2ND  RESIDUAL  URT  AIK)  IS  DA(K,I)  AT  THE  VERTICES 
C 11*1,2,3) 

S TOR  A*R*CP 1 
S T OR  B = R*T 

CAC1,  I )*STORA*X*STORB*Ut  OA ( 2, I »«STORA*Y*STORB*V*  DAI3,I)«STORA 
STOR*STOR6*UY 

CAI4,  IMSTORAX*  OAI  5,  I l*STOR*Y$  DA(6,I)=ST0R 

S T OR  A*CR  *T 

STOR*STORAARP 

CA(7,I)*ST0R*X4T$  DA  I 8, I ) »STOR*Y*  DA(9,I)*ST0R 

STOR*STORI*R  I 

CAI 10, I)  = STOR*X*  CAI 11, I )»STOR*Y$  OA 1 12 , I ) * S TOR 
C 2NC  DERIVATIVES  OF  D ( DA A I L , K , I ) ) ARE  NOT  TO  BE  COMPUTED  IF  INOI  = l 
IMINCI.FC.il  GO  TO  35 

C DAA(L.K.I)  IS  THE  CERIVATIVE  OF  THE  2ND  RESIDUAL  0 WRT  AIL)  ANO  AIK)  AT  THE 
C VFKTICES  (1*1,2,31 
RTT2=R*T  *T 
CRTT2*CR  *T  *T 

CAA( l, 1, I )»2.0*RTT2AX*  OA A ( 1 , 2 , I ) *R T T2*Y$  DA A ( 1 , 3 , 1 ) *RTT2 
CAAI 1,4, I)*0.0»  OAAI 1,5, I )*0.C*  DAA I l , 6, I ) =0 . 0 
ST0RA*ICR1*XaTAU)*T 
STOR=RP*STOR A 

CAAI 1,7, I )*STOR*X$  0 A A ( 1 , 8 , I I *STOR • Y«  DA A 1 1 , 9 , I ) *STOR 
STOR*RI ASTORA 

CAAI 1, 10, I )*STOR*X$  DAAI 1,11,1  )=STORAY»  DAA (1,12,1) *STOR 
CAA(2,2,1 )*0.0t  DAA(2,3,1)*C.C 

C A A ( 2 , A, I )*RTT2*Xt  DAAI2.5, I )*RTT2*Yt  DAA ( 2 , 6 , I ) *R TT2 

ST0RA*(CR1*Y»T*V)*T 

ST0R=RP*ST0RA 

CAA(2,7,I )»STOR*X*  DAA ( 2 , 8 , I ) «STOR*Y*  DAA ( 2 ,9 , I ) *STOR 
STOR*RI*STORA 

CAAI 2, 10, I )*STOR*X$  OAAI 2,11,  I )*STOR*Y$  DAA 1 2 , 12 , 1 ) -STOR 

CAAI 3, 3, I 1*0.01  OAAI 3,4, I )*0.0*  DAA I 3, 5, I I *0 .0$  DAA 1 3 ,6 , 1 ) *0.0 

ST0R*CR1TT*RP 

CAAI 3,7, I )*STORAX*  D A A I 3 , 8 , I) * STOR* Y*  OAA 1 3 , 9 , I I *STOR 
STOR*CRi T T*R  I 

CAAI 3,10, I l-STORAX*  0 A A I 3 , 1 1 , I)  * STOR * Y$  DAA I 3 , 12 , 1 I -STOP 
CAAI 4,*,  I 1*0.0*  D A A I 4 , 5 , I )*O.Ct  DA A I 4 , 6 , I I *0 . 0 
STORA=UYTTT*X 
S TOR*RP*STOR  A 

CAA(4,7,I )»STOR*X$  C A A 1 4 , 8 , 1 I * STOR ♦ Y t OAA 1 4 , 9 , 1) » S TOR 
STOR*RI*STORA 

CAAI 4, 10, I )*STOR*X$  OA A I 4 , 1 1 , I ) *STOR«Y«  OAA 1 4 , 12 , 1 ) -STOP 
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CAA( 5,5, I 1-0.0*  0AA<5,6, I 1-0.0 

ST0RA-UYTTT*Y 

ST0R-RP*ST0RA 

CAAI5.7, I l-ST0R*XS  0 A A < 5 , 8 , l I ■ S TOR* Y*  OAAt 5,9, I I -STOR 
ST0R-RI*ST0RA 

CAA(5, 10, I t«ST0R*X*  0AA(5,11,1)-ST0R*Y*  OAA  < 5 , 12  , I I -STOR 

C A A ( 6,6, I 1-0.0 

STOR-RP*UYTTT 

C A A ( 6 , 7, I |-ST0R*X$  0 AA ( 6 , 8 , I ) - ST0R*Y*  OAA (6 , 9 , 1) -STOR 
STOR-R I *UYTT  T 

CAA(6,lO, I )-STOR*XS  DAA(6, 11, I l«STOR*Y*  OAA ( 6 , 12 , 11 -STOR 
STOR A-CRT T2*X 
STOR-RPP*STOR A 

CAA< 7,7, II»STOR*X*  0 AA C 7 , 8 , I I « STOR*Y*  OAA ( 7 , 9 , 1) -STOR 
STOR-RPI*STORA 

C AA ( 7 , 1 0 , I )-ST0R*Xt  0AA(7,11,I)-ST0R*Y$  OAA ( 7 , 12  , I ) -STOR 

ST0RA-CRTT2*Y 

STOR-STO«A*RPP 

C A A ( 8 , 8, I l»ST0R*Y$  DA A (8, 9, I I -STOR 
ST0R-ST0RA*RP1 

C A A ( 8 , 1 0 , I )-STOR*X(  DAA( 8,11,  I l=STOR*Y$  OAA ( 8 , 12 , 1) -STOR 
C A A ( 9 , 9 , I l-RPP*CRTT2 
STOR-RP I *CRTT  2 

CAA { 9, 10, 1) =ST0R*X*  0 A A ( 9 , 1 1 , I I -ST0R*Y*  OAA ( 9 , 1 2 , 1 I -STOR 

ST0RA=CRTT2*RII 

STOR-X*STORA 

CAA( 10, 10, I l = STOR*X$  DAA ( 10 , l 1 , l I »STOR*Y*  OA A ( 1 0 , 12 , 1 I - STOR 
STOR=Y*STORA 

CAA( 11,11,  I l=STOR*Y*  OAA I 1 1, 12, 1 I -STOR 
CAAI 12, 12, I l-STORA 
CO  40  P=i,ll 
KP 1-Pa l 

CO  40  L-KP1.12 
40  C A A ( L , M, I )«CAA(P,L,I  I 
35  CONTINUE 

C FOR  X-POPENTUP  ECUATION  THE  CHARACTERISTIC  OUANT I T Y I CPCF I IS  VELOCITY 
C —TAKE  SOUNC  SPEEC  AT  T-0  ANO  CENTER  NODE 

CALL  SNOSPOIAPI  1),AE(1),DMCFI 
CO  300  IC-1,3 
Cl  IC  I-WTSI2  1*01  IC  l/OPCF 
CO  300  JC- 1 , PU A 

CA  (JC, ICI-WTSI 2I*CA(JC, ICI/DPCF 
CO  300  KC-l.PUA 

>00  CA A(KC, JC, IC  >-WTS<  2 l*DAA(KC, JC,  IC l/OPCF 
RETURN 

C RFC, IN  CALCULATIONS  FOR  K-3  (Y-POPENTUP  EON. 1 

45  VOX- I I AV!2I-AV( 31 l*AY( 1 1- ( A V ( 1 l-A V ( 3 I I *A Y I 2 I * ( A V ( 1 1 - AV 1 2 I »*AYI3I  I 
1 /CET 

VOY-(-( AVI  2 I -AVI  3 I »*AX(1 »♦( AVU l-AVI 31 l*AXI2 I- ( A V ( 1 1 - A V ( 2 I »*AX(3I I 
1 /CET 

POY-I-I AP(2)-AP(3H*AX(1I*(AP(1I-AP(3I»*AX(2I-(AP(1I-AP(2I I* AX  131  I 
l /CET 

VX«V0X*A(4I*T$  VY* VOY* A ( 5 »*T 

TTT=T*T 

CR5- 1 .0*  T *VY 

CR5TT-CR5*T 

CO  50  1*1,3 

X-AXI I I ( Y- AY ( I I 

R«  ARHO ( I , 1 1 1 RP-ARHO 11,21*  RI-ARH0(I,3I 


67 


RF  P- ARHO I I , 4 ) $ R I I -ARHOI 1,51$  RPI -ARHOI I . 6 I 
UT-A(1)4X«A(2)4Y«A(3) 

VT-A(4|AX4AI5)*Y+A<6> 

U-AUI  I IaTAUT 
V>AV(II«r«VT 
CR-VTaU*VXaVAVY 

C THE  THIRO  RESIDUAL  D AT  THE  VERTICES  OF  THE  TRIANGLE  (1-1,2,31  IS 
Cl  II-RACR*P0YaAI8|AT 

C THE  CERIVATIVE  OF  THE  3RD  RESIDUAL  *RT  A I K I IS  DA(K,I)  AT  THE  VERTICES 
C ( 1 = 1,2,31 

STOR-R*VX*T 

C A ( I , 1 1 -STOR*X$  DA( 2, I )-STOR*YS  DA(3,II=STOR 

ST0RA-RACR5 

STORB-RAT 

CA(<i,  I l>X*STORA«STORB*Ut  OA( 5, li -YASTORAaVASTORBI  DA  1 6 , 1 1 -STORA 

STORA-CRAT 

STOR»RP*STORA 

CA  I 7,  I l»STOR*X$  DA( 8, n-STORAV  + Tt  DAI9,II»STOR 
STOR-RI *STOR  A 

CAliO, I l-STORAXS  CAU1, 1)-STORAV$  DA  1 12 , I I - S TOR 
C 2NC  DERIVATIVES  CF  D (OAA(L,K,I>)  ARE  NOT  TO  BE  COHPITED  IF  INDI-1 
IF ( INC  I .EC. 1 ) GO  TO  50 

C C A A ( L , K , I ) IS  THE  CERIVATIVE  OF  THE  3RD  RESIDUAL  D WRT  AILI  AND  AIK)  AT  THE 
C VERTICES  (1-1,2,31 
RTT?-R*TTT 
STCRB-TTTAX 

CAA( 1,1, I 1-0.0$  DAA( 1,2, II-C.CS  OA A ( l , 3, I I >C . 0 
CAAI 1,4, I )»R*STORB$  C A A ( l , 5 , 11  -C . 0$  DAA ( 1 ,6 , I) >0 . 0 
STORA-VX ASTORB 
S TOR»RP*S  TOR  A 

CAAI 1,7, I l-STORAXS  0 A A ( l , 8 , 1) * S TOR A y $ DAA 1 1 , 9 , I I -STOR 
STOR-RI*STORA 

CAAI 1,10, ll»STOR*XS  DAAI 1,11, I l-STORAYS  DAA ( 1 , 12  , II -STOR 
STOR  B -TT  T *Y 

CAAI 2,2, I)-0. 0«  DAAI 2, 3,  I l-C.C 

CAAI 2, A, I >»R*STORB$  C A A I 2 , 5, I ) -C . 0*  OAA 1 2 , 6 , 1 ) -C .0 

STORA-VXASTORB 

ST0R-RP*ST0RA 

C A A I 2 , 7, I )-STOR*X$  DAAI 2, 8, I »«STOR*Y$  DAA I 2 , 9 , I) « STOR 

stor-ri*stora 

CAA(2«10«II»ST0R*X$  CAAI 2,11, I l-STOR*Y$  OA A I 2 , 1 2 , 1 ) -STOR 
CAAI 3, 3,  I 1-0.0 

CAAI 3,4, I I-RTT2I  CAAI 3,5,  I 1-0. 0$  DAA ( 3 , 6 , I I - C. 0 

STORA-VX-TTT 

STOR-RPASTORA 

CAAI 3,7, I l-STORAXI  D A A ( 3 , 8 , 1 I » S TOR a y $ D AA I 3 , 9 , 1) - STOR 
STOR  *R I *STOR A 

CAAI 3,10, D-STORAX*  C A A I 3 , 1 1 , I) -STOR AY$  OAA I 3 , 1 2 , 1 I -STOR 
CAAI 4,4,1  )-0 .Ot  DAA(4,5, I )-RTT2*X$  C AA I 4 , 6 , I ) -0. C 
ST0RA»XACR5TT«TTTau 
STOR-RPASTORA 

CAA(4,7, I )-STOR*XS  D A A I 4 , 8 , l I - S TOR A Y$  OAA 1 4 ,9 , 1) -STOR 
STOR-RI *STOR A 

CAAI4.10, D-STORAXI  DAAI 4, 11, I l-STORAYS  DAA ( 4 , 12 , 1) -STOR 
CAAI 5,5, I l-2.0*RTT2*Y«  DA A I 5 , 6 , 1 ) -RTT2 
ST0RA-YACR5TTATTTAV 
STOR-RPASTORA 

CAA(5,7, I l-STORAXS  DA A ( 5, 8, I I -STOR A Yi  OAA 1 5 ,9 , I) -ST  OR 
STOR  *K I AS  TOR  A 

CAAI 5, 10, I l-STORAXS  C A A I 5 , 1 1 , I I * STOR A Y*  OAA I 5 , 12 , 1 I -STOR 
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CAA(6,6,  I 1-0.0 
ST0R-CR5T  T *RP 

CAAI6.7, I l-ST0R*XS  0AA(6,8, I l«ST0R*Yt  OAA ( 6 . 9 , 1 I - STOR 
ST0R«CR5TT*R I 

CAA(6.10, I »-STOR*X*  OAA( 6, 1 1 . I I - STOR*Y»  OAA ( 6 . 12 . 1 1 -STOR 

STORB»CRM*T 

STORA-STORB*X 

STOR-RPP*STORA 

C A A ( 7 , 7, I l»STOR*X$  DAAI 7»8» I I -STORAYA  DAA I 7 ,9 , 1 I -STOR 
0AA(7,10»II-ST0R*X$  DA A I 7 , 1 1 » I I -STOR* Y*  OAA ( 7 , 1 2 , 1 I -STOR 
STORA-STORB*Y 
STOR-STORAARPP 

C A A ( 8 « 8, I )«  STOR*Y  $ CAA(8,9, I I -STOR 
STOR-RPI*STORA 

CAAI8.10.I l-STORAXJ  0 AA( a , 1 1 , 1 1 -STOR*Y$  OAA ( 8 . 12  • 1 1 -STOR 
CAA(9,9,I  | «RPP*STORB 
STOR- RP I AST ORB 

CAAC9, 10, I)*STOR*Xt  DAAI 9, 1 1 , I I -STORAYA  OAA ( 9 , 12 , 1) -STOR 

STORA-RI I*STORB 

STOR-STORAAX 

CAA< 10,10, II-STOR*X$  DAA I 10 , 1 1 , I) -STCRAVA  OAAI 10.12.1  I- STOR 
STOR-STOR  A ay 

CAA( 11.11. I l-STORAYA  DAA( 11,12, II -STOR 
CAA( 12,12,1 l-STORA 
CO  55  P-1,11 

KP1-P41 

CO  55  L-KP1.12 
55  CAA(L,P, I l«DAA(P,L, II 

50  CONTINUE 

C FOR  V-POPENTUP  ECUATICN  THE  CHAR AC  TER  I ST IC  OUANT I T Y I DPCF I IS  VELOCITY 
C —TAKE  SOUNC  SPEEC  AT  T-0  ANO  CENTER  NOOE 

CALL  SNOSPO(AP(l>,AE(ll,OPCFI 

co  aoo  ic-i.a 

C( 1C  l-WTSI 3I*D( IC  l/CPCF 
CO  400  JC-l.PUA 

CA  (JC,ICI-NTS(3IA0A(JC,ICI/DPCF 
CO  400  KC-l.PUA 

400  CAA(KC,JC,lCI=taTS(3IADAA(KC,JC, ici/dncf 
RETURN 

C MFC  IN  CALCULATIONS  FOR  K-4  (INTERNAL  ENERGY  EQN.I 

60  UOX-I  ( AU(2I-AU(  31  |AAY( ll-(AU(  1 l-AU  ( 3 I I *A  Y I 2 I ♦ ( AL  ( 1 1 -AU  (2  I I *AY  (31  1 
1 /CET 

V0Y=(-(AV(2I-AV(3I I*AX(IIa(AV(1I“AV(3I I *AX ( 2 I- ( A V ( 1 l-AV (2 l I ♦ AX (31  I 
1 /CET 

EOY=(-( AE ( 2 l-AE ( 3 I I A AX ( 1 1 4 ( AE ( l l-AE ( 3 I I A AX ( 2 I- ( AE ( 1 1- AE (2  I I AAX ( 3 I I 
1 /CET 

EOX-I  (AE(2  I - AE  ( 31  I • AY( 1 I - ( AE ( 1 l-AE ( 3 1 I * A V ( 2 I ♦ ( AE ( 1 1 - AE (2  I I * AY ( 3 II 
1 / CET 

EX-F0X*A(10I*T$  EY-E0V4A( 11IAT 

UXPVY-U0X«V0Y4T*( A( 1 |4A( 51  I 

TTT=T*T 

CO  65  1-1,3 

X- AX ( I I I Y* AY ( I I 

R»  ARHO( I , 1 I A RP -ARHO ( I , 2 1 1 RI»ARH0(I,3I 

RF  P = ARHO (1,411  RII-ARHOI 1,51*  RP I - ARHOI I , 6 I 

U-AU( I >♦?«( A( 1I*X4A(2I«Y4A( 31  I 

V- AV ( II*T*(A(4l*X*A(5)*Y4A(6t  I 

P-AP( 1 I«T*(A(7I*X4A(8)*Y4A(9I  I 

CR-A( 10 1 4X4 A ( 11 I*Y4A(12|4U*EX«V*EY 


69 


r 


L 


CRTT  *CR*T 

CR10«X*T*U 

£R11»Y*T*V 

C THE  FOURTH  RESIDUAL  D AT  THE  VERTICES  OF  THE  TRIANGLE  (1*1,2,31  IS 
C ( l)-R*CR*P*UXPVY 

C THE  DERIVATIVE  OF  THE  4TH  RESIDUAL  WRT  AIK)  IS  DA(K,I)  AT  THE  VERTICES 
C ( 1*1, 2, 3) 

STORB«P*T 
STOR  A*R*T 
STOR*$TOrtA*EX 

CAIl,  I >»STOR*X*STORB$  OA ( 2 , I ) «STOR*Y$  DA(3,II-ST0R 
STOP*STORA*EY 

CAI4,  I l=STOR*XI  DA (5, I )*STOR*Y*STORB$  DA|6,I)-ST0R 
STOR»RP*CRTT«UXPVY*T 

C A ( 7,  I ) *STOR*X$  OAI8,  I l-STOR*Y»  DA  I <) , I I *STOR 
STOR  = R I *CRTT 

DA (10, I )*ST0R*X*R*CR10t  DA ( 1 1 , I) »ST0R*Y*R*CR 11  * DA ( 12 , 1 ) *S TOR*R 
C 2ND  DERIVATIVES  OF  D ( C AA ( L ,K , I ) ) ARE  NOT  TO  BE  CCHPUTED  IF  INOI*l 
I F ( INDI.EC.l)  GO  TO  65 

C DAA(L,K,  I I IS  THE  DERIVATIVE  OF  THE  4TH  RESIDUAL  D hRT  AIL)  AND  AIK)  AT  THE 
C VERTICES  11*1,2,3) 

RTT2=R*TTT*  RTT 2TY»RTT2*Y$  RTT2TX*RTT2*X 
R PTT  = RP*T  $ R ITT-R l*T 
CR 1 l TT=CR 1 1 *T 
CRTT  2=CR  T T *T 

CR10RPT=CR1C*T*RP$  CR  1 1TRP-CR1 1 TT*RP t CRl 1R I T *CR 1 l TT*RI 

STORB*TTT*EX 

STORA=STORB*X 

STOR*TTT-»STORA*RP 

CAA(1,1,I)-0.0S  DAAI 1,2,  I 1-0. Ot  DAA ( l , 3, I I =C .0 
CAAI 1 , A,  I 1*0.0$  DAAI 1, 5, I l-C.CS  DAA (1,6, I ) =0 . 0 
CAAI 1,7,1 >*STOR*X»  DA A 1 1 , 8 , 1 I * STOR* Y$  DAA 1 1 , 9 , II - STOR 
STCR  = STORA*R  I 

CAAI 1, 10, 1 )*RTT2TX*ST0R*X$  OA A ( 1 , 1 1 , I I *STOR* Y$  DA A 1 1 , 12 , I ) *STOR 
CAAI 2,2,  I 1*0.0$  DAAI 2,3, I >*C.C 

CAAI2.A, I 1*0.0$  DAAI2.5, I l-O.Ct  DAA ( 2 , 6, I I *0 .0 

STOR A*Y*ST0RB 

STOR=RP*STORA 

CAA(2,7,I  1»ST0R*X8  DAAI2,8, 1 1«ST0R*Y$  DAA 1 2 , 9 , I ) «STOR 
STOR*RI*STORA 

CAAI2.10, l )«ST0R*X*RTT2TY$  DA A I 2 , 1 1 , I I *STOR* Y $ D A A 1 2 , 12 , I I = STOR 
CAAI3, 3,  11*0.0 

CAAI 3,4, I 1*0. 0$  DAAI 3,5, I l-C.C*  OAA I 3, 6, I ) *0 . 0 
STCR=STORB*RP 

CAAI 3,7, I l-STOR*X*  0 A A ( 3 , 8 , 1 I * S T0R*Y«  DAA <3 , 9 , 1 I * S TOP 
STOR»STORB*R  I 

CAAI 3, 10, I )*ST0R*X+RTT2t  OAA ( 3 , 11 , I ) *STOR*Y t OAA  I 3 , 12 , 1 I * STOR 

CAAI A, A, I 1*0. 0$  DAAI A, 5, I ) *0.C$  DAA I A, 6 , I ) *C  .0 

STORB*TTT  *EV 

ST0RA*ST3RB*X 

STOR*STORA*RP 

C A A I A , 7, I |*STOR*Xt  0AA(4,8,  I l«STOR*Y*  DAA < 4 , 9 , II *STOR 
STOR*STOR A*R  I 

CAAI 4,10, I )*STOR*XI  CAAI 4,11, I l*ST0R*Y*RTT2TX$  0 A A (4 ,12 , 1 I *STOR 
CAA(5,5, I 1*0.0$  DAAI 5,6, I l-C.C 
STOP A=STORe*Y 
STOR*STORA*RP*TTT 

CAAI 5,7, I l*ST0R*Xt  DA A ( 5 , 8 , 1 I -STOR* Y$  DAA 1 5 , 9 , I I * STOR 
STCR*STORA*R I 

CAAI 5, 10, I l*STOR*Xf  DAA ( 5 , l l , I ) -RTT2 TY*STOR*Y»  DA A 1 5 , 12 , 1 I -STOR 
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CAA(6,6,  I 1*0.0 
STOR*STORB*RP 

CAA(6,7, I >-STOR*X$  OAA( 6, 8, I) «STOR*Y « DAA I 6 ,9 , I I -S TOR 
ST0R*ST0Rfi*R  I 

CAA(6,10»I)*ST0R*X$  DAA ( 6, 1 1 , I I «RTT2*STOR* Vt  DA A ( 6 , 12 , 1 I- S TOR 
ST0RA-CRTT2AX 
ST0R*ST0RA*RPP 

C A A ( 7 , 7, I l*STOR*X$  DAA ( 7 , 0 , I)  »STOR*Y*  OAA ( 7,9,1  I »ST0R 
ST0R=ST0«A*RP  I 

CAA( 7, 10,  I l*(STOR«CR10RPT)*X$  OAA ( 7 , 11 , I ) »ST0R*Y ♦CR11TRP*X 
CAA( 7,12, I )*STOR  + RPTT*X 
STOR A»CRT  T2* Y 
ST0R*ST0RA*RPP 

C A A ( 8 , 8,  I I * S TOR  *Y  $ OAA (8,9, I l*STOR 
STOR*STORAARP  I 

CAA( 8,10, I )«ST0R*X*CR10RPT*Y$  0 A A ( 8 , 1 1 , II  * ( S TOR + CR UTRP ) * V 
C A A ( 8 , 12 , 1 1 »STOR-»RPTT*V 
C A A ( 9, 9, I )-CRTT2*RPP 
ST0R=CRTT2*RP I 

CAAI9.10, I )*ST0R*X>CR10RPT»  DA A ( 9, 1 1 , I I »ST0R* Y*CR1 1 TRP 
CAA(9,12,  I l = STOR*RPTT 
ST0RA=RII*CRTT2 
ST0R*ST0RA*X4RITTACR10 

CAA( 10,10, I > = ( STOR^RITT*CR10l*Xt  OAA ( 10, 11 , I I * S TCR* Y+CR 1 1R I TAX 
C A A ( 10,12,1 l-STOR«RITT*X 
STOR»STORA*Y 

CAA( 11,11, I I=ST0R*Y*2.0*CR11RI T*Y$DA A ( 11,12,1 I»ST0R*CR11RIT*RITT*Y 
CAA( 12, 12,11 =ST0RA*2.0*R ITT 
CO  70  P*l,ll 
KPi=M*l 

CO  70  L*KP  1 * 12 
70  CAAIL.M, I )«DAA(P,L,l  I 
65  CONTINUE 

FOR  ENERGY  EQUATION  TEE  CHARACTERISTIC  QUANTI TY ( CPCF I IS  INTERNAL  ENERGY 
AT  T*0  ANC  CENTER  NOCE 
CPCF= AE ( 1 I 
CO  500  I C*  1, 3 
C( IC I«WTS(A1*D( IC l/CPCF 
CO  500  JC*1,PUA 

CA  ( JC, IC  >*WTS( A l*DA( JC, IC  l/DPCF 
CO  500  KC*1,PUA 

500  CAA(KC,JC, ICI-WTS(A)*DAA(KC,JC, ICI/D MCF 
RETURN 

101  PRINT  98 , PUA 
STOP 

98  FORMAT ( IH, • PUA» • , 1 5 , • BUT  SLBROLTINE  OSUB  IS  WRITTEN  FOR  PUA=12,ST 
1CP ' I 

99  FORMAT  ( 1H,»  K*', 15, ‘BUT  SUBROUTINE  DSUB  IS  WRITTEN  FOR  K>1,2,3  0 
1R  A.  STOP • I 

ENC 

SUBROUTINE  AREQST(NUMOND,NOEQ,PUA,NDPPRP,GAUST,X,Y,VALS,B,RHOI 
C COMPUTES  CENSITY  ANO  ITS  DERIVATIVES  AT  ALL  NODES  OF  A PRISM 
C AT  A GIVEN  GAUSSIAN  TINE 

Cl  MENS  ION  X ( NUPONC  I , Y ( NUMCNO I , V AL  S ( NOF Q ,NUPOND  I ,8 (PL A I ,RHO ( NUPOND, 
110),R(10I 
CO  50  K*i,NCPPRP 

P<VALS( 3,K l«GAUST*(B(7l*X(K l«B(8l*Y(K|AB(9) I 
E*V*LS(A,K)aGAUST*(B(10)*X(K»*B(11I*Y(K>*B(12)I 
CALL  EQNST ( P , E ,R  I 
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CO  50  1=1,10 
50  RHCIK.I l=R<  I 1 

C IF  NUMBER  OF  TRIANGLES  IS  LESS  THAN  MAX  NUMBER  ALLOWED  IN  A PRISM,  SET 
C RFMAINCER  OF  RHO  TO  ZERO 

IFINDPPRM.EC.NUMONO)  GO  TO  101 
J=NCPPRM«1 
CO  60  K = J ,NUMONO 
CO  60  1=1,10 
60  RHO ( K , I (=0.0 

101  CONTINUE 
RETURN 
ENC 

SUBROUTINE  ECNST(P,E,R) 

C COMPUTES  CENSITY  FOR  A NOBEL  ABEL  GAS  AND  ITS  FIRST  9 PARTIAL  DERIVATIVES 
C AT  A GIVEN  GAUSSIAN  TIME 

C HAS  TWO  PARAMETERS  GAMMA  AND  ETA  (UNITS  FOR  ETA  ARE  M**3/KG) 

C I MENS  ION  R ( 10  I 

DATA  GAMMA, ETA/1. 40, 0.0/ 

G = G AMM A-  1.0 
C=G*F*FT  A*P 
IFIC.EQ.O.OICOTO  101 
CO=C*C 
DCC=D*D*C 
CCCC=CD*CC 
C hi  1)  = CENSITY 
Rill =P/C 

C M2)=FIRST  DERIVATIVE  OF  Rill  WRT  PRESSURE 
R l 2 I =G*E/CC 

C *m  = Fl«ST  CERIVATIVE  OF  R(ll  WRT  INTERNAL  ENERGY 
R(3I=-G*P/CC 

C R(4l«CERIVATIVE  OF  RI2I  WRT  PRESSURE 
RI4I»-2.0*ETA*G*E/C00 

C R ( 5 I = DER  I VAT  I VE  OF  R(3)  WRT  INTERNAL  ENERGY 
P(5)=2.0*G*G*P/CDC 

C R«6l*CERIVATIVE  OF  R(2)  WRT  INTERNAL  ENERGY  OR  OF  R(3)  WRT  PRESSURE 
R(61=-G*(G*E-ETA*P1/CC0 

C R ( 7 1 IS  THE  3RC  PARTIAL  OF  OENSITY  WRT  PRESSURE 
R I 71=6.0* ETA* ETA*G*E/CCDD 

C R(BI  IS  THE  3RD  PARTIAL  OF  DENSITY  WRT  INT.  ENERGY 
R(8I=-6.0«G*G*C*P/CCCD 

C R(4)  IS  THE  3RD  PARTIAL  OF  DENSITY  WRT  INT.  ENERGY  AND  PRESSURE  SQUAREC 
R(9I*2.0*ETA*G*(2.0*G*E-P*ETAI/DDOD 

C RIIOI  IS  THE  SRC  PARTIAL  OF  OENSITY  WRT  PRESSURE  AND  INT.  ENERGY  SQUARED 
R( 10  l=2.0*G*G*(G*E-2.0*ETA*P  l/DCOD 
RETURN 

lOi.  PRINT  102 

102  FORMAT! 1H,//,2X, 'EQUATION  OF  STATE ,OENSI TY I PRE S SURE , I NTERN AL  ENERG 
1Y1.IS  WRONG  SINCE  DENOMINATOR  IS  ZERO.  STOP  SUBR.  EQNST'I 

STOP 

ENC 

SUBROUTINE  SNDSPOI P, E , SS  I 

C COMPUTES  SOUNC  SPEED  FOR  A NOBEL-ABEL  GAS  WITH  PARAMETERS  GAMMA ,ET A 
C UMTS  FOR  FT  A ARE  M**3/KG 

CATA  GAMMA, ETA/1. 40, 0.0/ 

G-GAMMA-1 .0 

S = GAMMA*I E*G*ETA*P  l**2/( E*G  I 

SS=SQRTF I S I 

RETURN 
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SUEROUT l NE  AREAIN(A,0,OET, VALLE) 

C COMPUTES  INTEGRAL  OVER  A TRIANGLE  OF  A PRODUCT  OF  TWO  FIRST  DEGREE  POLYNOMIALS 
CIMENSION  A( 3). 01 3) 

C * A I 1)*0(  1)+A(2)*B(2)+A( 31*0(3) 

C“  A I 1)*B(2)+A(2)*B( 1 ) 

E*A(l)*Bl3)«A(3)*0il) 

E=A(2I*B<  3 1 -» A I 3 ) * 0 I 2 » 

VALUE*ABSF<  CET)*(C/12.0*(C*E*F)/25.0) 

RETURN 

ENC 

FUNCTION  CETfcRMlX.Y) 

C COMPUTES  I 1 , 1 , 1 ) COT  I Y CROSS  X) 

CIMENSION  X<  3 ) • Y ( 3 ) 

CETERM=X(1)*IY(3)~Y(2))*Y(1)*CX(2)-X(3))*(X(3)*Y(2)_X(2)*Y(3II 

RETURN 

ENC 

SUBROUTINE  INI TZR I NOEQ , MUA , NONUM , A ) 

C REPLACES  INITIAL  CATA  USED  AT  T IME=OLO  TIME  WITH  INITIAL  DATA  TO  BE  USED  AT 
C T I ME  = OLC  T I MF  *CT 

COMMON/BLOK3/VARTI A, 2000) 

CIMENSION  A ( MU A ) 

C ASSUMPTION  PARAMETERS  All)  ARE  EQUALLY  SPACED  AMONG  UNKNOWNS 
MC IV*MUA/NOEQ 
CO  20  K=i,NCEQ 

C SAVE  ONLY  LAST  PARAMETER  IN  EXPRESSION  FOR  EACH  FLOW  VARIABLE — INITIAL  VALUE 
C OF  THE  FLOW  VARIABLE  AT  NEXT  TIME  STEP 
20  VART IK, NONUM  )*A(MC  IV*K  ) 

RETURN 

ENC 

SUBROUTINE  V ALNEW ( NUMOND , NOEQ , MUA , KNEW , DT , NONUM  , VALS  , A ) 

C SUBSTITUTE  NEW  VALUES  OF  THE  FLOW  VARIABLES  AT  T I ME«CLD  TIME+DT  INTO  MAIN  INFO 
C MATRIX  VAR  AT  LEVEL  KNEW 

COMMON/0LOK2/VARI 2, A, 2000) 

CIMENSION  VALS (NOEQ, NUMOND ) , A I MUA) 

C ASSUMPTION  PARAMETERS  All)  ARE  EQUALLY  SPACED  AMONG  UNKNOWNS 
MC I V =MU A/NOE  C 
CO  20  K=l,NOcQ 

20  VAR(KNEW,K,NCNUMl=VALS(K,  1 ) *0 T* A I MD I V*K ) 

RETURN 
ENC 

SUBROUTINE  OLTN06  I NUMBER , T I ME , KNEW ) 

UUTN06  IS  AN  OUTPUT  ROUTINE  WHICH  PRINTS  THE  KNEW  LEVEL  OF  VAR  AND  THE 
CENS  ITY 

COMMON  /PL0K2/VAR<2,5,2000) 

PRINT  500 

<.00  FORMAT  I 1F,//,20X,  ‘MATRIX  VAR,  AND  DENSITY*) 

PRINT  NO l 

501  FORMAT  I1F,/,5X,*N0DE*,9X,*X  VEL*.15X,*Y  VEL 1 5 X ,* PRESSURE ', 15X , 

1 • INT  ENG*, 11X, 'DENSITY*  ) 

CO  15  J*  l , NUMBER 

CENS*2.5«VAR(KNEW, 3, J ) /VAR  I KNEW, 4, J ) 

15  PRINT  502,J,VAR(KNEW,1,J),VAR(KNEW,2,J),  V AR  ( KNE  W ,3  , J ) , 

1 VARIKNEW.A, J),DENS 

502  FORMAT! IH,3X, l 5,51 5X.1PE 16.8)  ) 
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RETURN 

ENC 

SUBROUT I NE  GRAPHIT  ( T , KNEW .NUMBER , INDEX, ORGT I M , BEGPOS  I 
C SUBROUTINE  GRAPHIT  TAKES  THE  VALUES  FROM  VAR  I KNEW , I , J ) NORMALIZES  THEM 
C PY  THEIR  CORRESPONCING  SHOCK  VALUES,  PRINTS  THEIR  VALUES,  AND  CALLS 
C SUBROUTINE  PLOTTO  TO  PLOT  THEIR  GRAPHS 
C THE  CURRENT  GRAPHIT  IS  WRITTEN  FCR  NODPOS=l,  NOTD-4 
C IMENS ION  RI200I,  CATP(4,2CC) 

COMMON/ BLOK  1/XY  ( 2, 20C0  ),  NR  ( <3 , 2000  )/BL0K2/VAR(  2,4, 20001 
COMPUTE  SHOCK  VALUES 

SUBSCRIPT  I CENOTES  QUIESCENT  STATE 
FORMULAS  ARE  FCR  GAMMA=1.4 
SHKSTH=  5.0 
OT IME*ORGT  IM 
OPOS=BEGPCS 
V1=0.0 

P1=1.01325E5 
Rl  = l. 225570786 
T1=288.15 

A1  = SCRT( 1 .4AP1/R1  ) 

SHVV  = Al*SCRT ( ( 1.0+6.0*SHKSTH)/7.0) 

SHKPOS=OPCS*SHVV*!T-OTlME  ) 

SHKV=5.Q*Al*( SHKST  H— 1.0)/SQRT(42.0A  SHK  STH  + 7 . 0 ) 

SHKR=R1»( l.C  + 6.0*SHKSTH)/( SHKSTH  + 6.0  ) 

SHKP=P1ASFKSTH 
SHKE=2.5*SHKP/SHKR 
N0CP0S=1*  1=0 
N0TC=4 
MCR=N0TC/2 
C RIGHTISH  POSITION 

51  IFINICR. GT. NUMBER IGOTO  34 
I = I ♦ 1 

R ( I ) =XYI 2 ,N I CR  I 

C A TP  I 1, l ) = VAR(KNEW, 3.NICR  I/SHKP 
C AT  P (2, I ) = VARIKNEW,2,NICR)/SHKV 
CATP  (4, I )»VAR(KN£W,4,NICR  l/SHKE 
CATP  ( 3, I I*  DATPI I, I )/DATP(4,I I 
C CENTER  POSITION 

MCC  = NICR  + N0TD  + 2 
IFINICC.GT. NUMBER  IGOTO  34 
1=1  + 1 

R I I»=XY(2,NICC) 

CATP  (1, I )*VAR(KNEW,3»NICC)/SFKP 
CATP  12, I l=VAR(KNEW,2,NICC)/SFKV 
CATP  (4, I )*VAR(KNEW,4,NICC)/SHKE 
CATP  13, I)=  DATPI 1, I I/DATPI4,I ) 

NICR=NICC+N0TD+1 
GOTO  51 

34  PRINT  198, T 

198  FORMAT! 1H , / / / , 10X, • T I ME  = • , 1PE16.8.10X, • PRE SSLRE , VELOC I T Y , DENS  I TV  , 
IENERGY  ARE  NORMALIZED  WRT  THEIR  SHOCK  VALUES') 

PRINT  197 

19/  FORMAT! 1H,//,4X,*  Y POS I T I ON • , 7X , • PRE SSURE  ' , 

19X , • Y VELOCITY', 8X, 'DENSITY', 9X,'INT.  ENG. • I 
PRINT  199, (RIKI, (CATP  ( J , K ) , J* 1 ,4 » ,K  = 1 , I ) 

199  FORMAT! 1H,(5! 3X , 1PE 16 . 8 I , / I I 
CALLPLOTTC(T»l»OATP,R,SHK PCS, SHKP,SHKV»SHKR, INDEX) 

RETURN 

ENC 
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SUBROUTINE  PLOTTO  I T ,NODP , DATAP , X, SHKPOS.PV , V V ,R V , INDEX ) 

C PLOTTO  USES  TEE  COHPLOT  ROUTINES  TO  PLOT  INFO  IN  OATAP 
C SEE  TECHNICAL  REPORT  ARDC  TR6  FOR  DEFINITIONS  OF  VARIABLES  ANO  SUBROUT 

C INES 

C I PENSION  SCF I A ) * TTL  8 1 10) 

C I PENS  ION  TTTLBUOI 
C I PENS  ION  TLBBI 10) 

CIPFNSION  SXI20),SYI20I 

Cl  PENS  ION  BUFFI  3000),  DAT  API  A,  200),  XI  200 1*  TLB  (10)  , SLA  I A)  .ARRI200I 
CATA  ISLAI  I >,  1*1,4), /'  PRESSURE)**  • Y VEL>*. 

I  * DENSITY)*,*  I ENERGY)'  / 

SCFI I ) = PV 
SCF I 2 ) = VV 
SCFI 3 ) =R V 

SCFIA)*2.5*SCFI1)/SCFI3) 

XPIN=2.8 
XP  AX*  3 . 3 
YP  IN=-0. 6 
YPAX=1.6 
X0=12.5 
VC*U.0 
XS*0.0A 
YS*0 . 2 

xe*io.o 

YB=-11.0 
X0R=2.fl 
V0R  = -0 .6 
CX=0.01 
CY*0.1 
H T =0 . 2 
XFAC*1.0 
YFAC*1 .0 
XPAGE=27.C 

PRINT  199,  XPIN,XPAX,YMIN,YPAX,XD,YO,XS,YS,HT,XCR,YOR,XPAGE 
. 99  FORPAT  I 1E,3X, *XMIN*» ,E15.8,3X, • XMAX« • ,E 15. 8 , 3X , *YNIN=* ,E15.8,3X, 

1 • YPAX=*,E15.B,/, 3X, *XD»* ,E15.8,3X, • YQ»  • ,E 1 5 . 8 , 3 X , • XS« • ,E15.8,3X, 

2 *YS**,E15.6,/  ,3X, *HT«* ,E 15 . B ,3X , • XOR* « ,E1 5. 8 , 3X , • YOR** .E15.8, 

3 3X, *XPAOE»* ,E15.8  ) 

SXI1)*SHKP0S 

SX  I 2 ) *SHK  POS 
SX I 3 ) *SHKP0S 
S X I A I *SHKP0S 
S Y I 1 ) * 1 • A 
SYI2)*1.0 
S Y I 3 ) *0 . 2 
S Y I A ) *-0  • A 
P*2 

CALL  PLTCCBI XPAGfc.M,  BUFF  I 1 ) , BUFF  I 3C00 ) ) 

CO  22A  JKL-l.A 
CO  50  JPN* 1 , NODP 
50  ARP  I JPNI =CATAPI JKL , JHN  ) 

YB*YB*18.C 

CALL  PLTCCS  I XB , YB , XOR , YOR , XS , YS ) 

PP*A 

CALL  PLTCCA  I OX , DY , XN IN, XPAX , YN IN, YP AX* MM ) 

S INA*  1.0 
C0SA*0.0 
XY*2. 72 
YX*0  • 2 
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CALL  PLTCCT(ET,SLA(JKU,  S INA , COSA , X Y , YX  ) 


C INDICATE  TINE  VALUE 
ENC00EI50,99,TLBI  T 
99  FORMAT  (•  TIME-  ',E15.8,1H>> 

SINA-0.0 
C0SA-1.0 
XY-2.76 
YX  — 0.95 

CALL  PLTCCT  < HT , TL B ( 1 > , S I N A , CO S A , X Y , VX ) 

C INDICATE  TEE  NORMALIZING  FACTOR 

SPSC-SCFI JKL  I 
ENC0DEI50,lIl,TTLei  SPSQ 

111  FORMAT  (•  NORMALIZING  FACTOR  ••,E15.8,1H>  ) 

YX— 1.15 

CALL  PLTCCTIET.TTLBI l > , S I N A ,CCS A , XY , YX » 

C INDICATE  TEE  TIME  LEVEL 

ENCODE  I 50, IC5»  TT  TLB  I INDEX 
105  FORMAT  I • TIME  LEVEL- ',15,  1EH>  I 
YX— O.B5 

CALL  PLTCCT 1 ET, TTTLB ( 11 , S IN A , CO SA , X Y , YX I 
ENCODE  <5C,112,TLBB) 

112  FORMAT  I • AESCISSA  IS  Y-COORCINATE  IN  METRES 
YX— 1.05 

CALL  PLTCCT I ET.TLBB! 1 ) , S I N A .CCS A , XY , YX ) 

CALL  LABEL A I CX,CY, XM I N , XMAX , YM I N , YMAX , XF AC , YFAC ) 

MNM-1  JNS-0  * IC-0 

CALL  PLTCCCIMMM.NS.XI  1),ARR( 1)  ,NOOP,IC,  XM I N , XM AX , YMI N , 

1 YMAX) 

NS-1U*  M=2  I N-5  * IC-0 

CALL  PLTCCC  (M.NS.SXU)  ,SY<1)  ,N,  IC,  XMI N , XM  AX  , YM  IN  , YM  AX  ) 

225  CONTINUE 

CALL  PLTCCP 

RETURN 

ENC 

* COMPILE  CISC, LABELA, ALL 


LIST  OF  SYMBOLS 


a 

a. 

1 

e(x,y,t) 

P(x,y,t) 

r 


t 

t 

At 


uCx,y,t) 

v(x,y,t) 

vr(x,y,t) 

x 

X 

y 

y 

Dk(x,y,t;a) 

E(a) 


F^a) 


vector  of  unknown  parameters 

. th  - -> 

1 component  of  vector  a 

internal  energy  per  unit  mass  [J/kg] 

pressure  [Pa] 

polar  radial  coordinate  [m] 

radii  of  annular  region  for  blast  wave  calcula- 
tion [m] 

time  [s] 

given  value  of  time  [s] 
time  increment  [s] 

Gaussian  quadrature  points  used  in  time  integra- 
tion [s] 

velocity  component  in  x direction  [m/s] 

velocity  component  in  y direction  [m/s] 

velocity  component  in  the  radial  direction  [m/s] 

Cartesian  spatial  coordinate  [m] 

given  value  of  the  x coordinate  [m] 

Cartesian  spatial  coordinate  [m] 

given  value  of  the  y coordinate  [m] 

nondimensional  residual  error  from  the  k^  differ- 
ential equation  at  a point  (x,y,t) 

total  residual  least  squares  error  over  a finite 
element 

the  first  partial  derivative  of  E(a)  with  respect 

to  a. 

1 


NOl:Q 

NUMTRI 

V 

P(x,y,t) 

u(x,y,t) 


number  of  differential  equations  to  be  solved 
simultaneously 

number  of  prisms  composing  a finite  element 
volume  of  a finite  element  [m3] 
density  [kg/m3] 
generic  flow  variable 
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corresponds  to  known  value  at  a given  time 
corresponds  to  value  at  the  shock  front 
corresponds  to  value  in  the  pre-shock  state 
corresponds  to  value  in  the  post-shock  state 
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